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Determination of Design and Operation Parameters 
for Upper Atmospheric Research Instrumentation 
to Yield Optimum Resolution with Deconvolution 

The subject NASA grant number NAG 1-804 was instituted 1 Aug 1987, originally 
for a three-year period, through 31 July 1990. No-cost extensions were granted so that the 
grant finally terminated on 31 Dec 1991. During the period in which the Principal 
Investigators were partially supported by the grant, Graduate Research Assistants and a 
Graduate Research Associate, Mr. Abolfazl M. Amini, were also supported on an 
intermittent basis. A good deal of grant-related research was performed by graduate students 
who were not directly supported by the grant, but were employed by the University of New 
Orleans Department of Physics as Teaching Assistants, or otherwise employed. 

In addition to the support provided by NASA Langley Research Center, two other 
sources of support were associated with this grant. Two funding increments through this 
NASA grant were provided by the U. S. Army Cold Regions Research and Engineering 
Laboratory (CRREL) for research performed by the Principal Investigators and an additional 
Investigator, Dr. Clyde Bergeron of the Department of Physics. A funding increment was 
also provided by NASA Marshall Space Flight Center to begin research into the analysis and 
prediction of the tethered satellite tether skiprope mode, which could occur in the 
NASA/Italy ASI TSS-1 experiment. 

The research for die basic NASA Langley grant and for the additional increments of 
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funding has been described in journal and proceedings articles, published abstracts, student 
Master’s theses, and reports, all of which are included in this final report document. We 
include both papers directly supported by NASA and the related research directed by the 
Principal Investigators which did not receive direct NASA support. A list of the journal and 
proceedings articles, research papers with published abstracts, and theses directly related to 
the subject of the initial grant is included in Appendix 1 , along with copies of the articles and 
abstracts. Appendix 1 is bound in this volume. 

Appendices 2, 3, and 4 contain the three student Master’s theses listed in Appendix 
1 . Each of these is bound separately and enclosed as a volume with this report. The two 
reports written to describe the research supported by CRREL are included as Appendices 5 
and 6, which are bound in this volume. Finally, the Engineering Notebook, Appendix 7, 
written for and submitted to NASA Marshall Space Flight Center, is also part of this report, 
and is included as an enclosed separate volume. 

Among all the research results reported in the Appendices, note should be made of 
the specific investigation of the determination of design and operation parameters for upper 
atmospheric research instrumentation to yield optimum resolution with deconvolution. As 
reported by G. Ioup et al (1988, 1989), a methodology has been developed to determine 
design and operation parameters for error minimization when deconvolution is included in 
data analysis. An error surface is plotted versus the signal-to-noise ratio (SNR) and all 
parameters of interest. Ins trumental characteristics will determine a curve in this space. The 
SNR and parameter values which give the projection from the curve to the surface, 
corresponding to the smallest value for the error, are the optimum values. These values are 
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constrained by the curve and so will not necessarily correspond to an absolute minimum in 
the error surface. 

During the period of this grant, the Investigators and their students have maintained 
frequent contact with the original technical monitor, Dr. George M. Wood, and the new 
monitor, Dr. Billy T. Upchurch, both of NASA Langley Research Center. This interaction 
has been immensely rewarding for both the Investigators and their students. We are very 
grateful to NASA Langley Research Center, not only for the funding of the research, but 
also for the interaction and research opportunities which have been provided to the 
Department of Physics at the University of New Orleans. We anticipate with pleasure 
continued association with the research staff of NASA Langley Research Center. 
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Autocorrelation estimation using constrained iterative 
spectral deconvolution 


Murali Ramaswamy* and George E. loupj 


ABSTRAC1 

Computing an autocorrelation conventionally pro- 
duces a biased estimate, especially for a short data se- 
quence. Windowing the autocorr Nation can remove the 
bias but at the expense of violating the nonnegativity of 
the corresponding power spectrum. Constrained itera- 
tive deconvolution provides a basis for improving an 
autocorrelation estimate by reducing the bias while 
guaranteeing nonnegativc definite ness. 

The length of the autocorrelation is increased in order 
to satisfy the nonnegativity constraints on the power 
spectral estimate. The constraints can also have signifi- 
cant effects on small, poorly determined values of the 
autocorrelation. The technique ;s applied to synthetic 
and real examples to show the improvements in the 
autocorrelation and power spectium which are possible. 

The method is reasonably sti ble in the presence of 
noise and it approximately preserves the area of the 
power spectrum. Comparison to (he maximum entropy 
technique shows that the iterative method gives power 
spectral resolution which is sometimes better and some- 
times not as good, but that ther : are cases for which it 
is the more desirable approach 


INTRODUCTION 

Autocorrelation and power spectral estimation are difficult 
when the segment of data available s of inadequate length. To 
reduce the variance in the power spectral estimate, auto- 
correlation windows (Blackman and Tukey, 1958; Geckinli 
and Yavuz, 1983) or some forms of spectrum averaging are 
used (Daniell, 1946; Bartlett, 1948; Welch, 1967). Bias in the 
autocorrelation due to missing lagged products, with more 
products missing as the lag increases (Cooley et al., 1970), is 
especially acute for short data segments, since the suggestion 
of Blackman and Tukey (1958) to use only 10 percent to 20 


percent of available lags is no longer practical. Short data 
sequences are commonly selected in designing filters for non- 
stationary or shift-variant data so that there will be approxi- 
mate stationarity or shift invariance within the window. This 
approach is often taken in processing seismic data. 

Maximum entropy (MEM) spectral estimation (Burg, 1975; 
Kanasewich, 1973) is the technique generally used for obtain- 
ing an improved autocorrelation and power spectrum for 
shop data sequences (Jurkevics and Wiggins, 1984). Maximum 
entropy spectral estimation assumes an autoregressive model 
for the data. When the method works, it w'orks very w’ell. 
However, it is sensitive to noise and is unstable for some data 
type:, (Lacoss, 1971; Chen and Stegan, 1974). It also does not 
return a true magnitude for the power spectrum (Lacoss, 1971 ; 
Johnson and Anderson, 1978). 

In this paper we suggest for improving estimates of auto- 
correlations and power spectra an alternative approach, which 
is not very sensitive to noise and assumes only that the data 
sequence is a truncated portion of a larger data set. 

The method begins with a biased estimate. It removes the 
bias using constrained iterative deconvolution of the power 
specirum. If run to convergence, the technique removes almost 
all the bias. For difficult data however, the iterations may be 
terminated before all bias is removed. The constraint guaran- 
tees the nonnegativc definite property of the autocorrelation 
by keeping the power spectrum nonnegative. The principal 
effect of the constraint is to extend the autocorrelation beyond 
the number of lags originally possible, but it can also improve 
the estimate in the original lag domain and reduce the effects 
of variance in estimating small autocorrelation values (Yoer- 
ger, 1978; Yoerger and Ioup. 1983). The constraints will most 
affect power spectra which have peaks with valleys near the 
zero baseline and least affect flat power spectra which rarely 
appmach the zero baseline. For all data, however, the pro- 
posed approach (Ramaswamy, 1985) removes the bias in the 
autocorrelation estimate and sharpens the power spectral esti- 
mate; i.e., the bias in the autocorrelation results from 
deemphasizing the larger autocorrelation lag values, which in 
turn smooths the spectral estimate. 
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We briefly discuss basic methods for autocorrelation and 
power spectrum calculations for truncated data sequences to 
show the bias which is inherent in the resulting estimate 
(Cooley et al., 1970). We then summarize iterative deconvolu- 
tion as it is currently practiced in the function domain, show 
how the technique is applied in the transform domain to im- 
prove the autocorrelation estimate, discuss practical consider- 
ations, and give some examples. 

AUTOCORRELATION ESTIMATION 

The autocorrelation of a time sequence x(t) that is N sam- 
ples long is often defined as (Cooley et al. f 1970) 

N- 1 

b(t) = 1/N £ x(t')x(t' - t). (!) 

t' = o 

The above formulation leads to a biased estimate of the auto- 
correlation because the number of lagged products going into 
the estimate of b(t) is a function of t. To correct this effect, the 
equation is rewritten for an unbiased estimate u{t\ 

N 1 

U(t) = 1/(N - U I ) £ x(Ox(r' - t). (2) 

V- 0 

Although the above equation results in reduced bias, this form 
can lead to a violation of the nonnegative definite property of 
autocorrelations of real sequences, corresponding to the non- 
negativity of the power spectrum. 

The biased estimate b(t ) is related to the unbiased sequence 
u(t) by the relation 

b(t) = u(t)w{t\ (3) 

where w(f) is a triangular or Bartlett window of unit peak 
amplitude at the origin (Cooley et al., 1970). The window can 
be expressed as 

w<0 = (N - 1 1 1 )/N, - N < t < N. (4) 

Uppercase notation is used to represent the transform vari- 
ables as follows: 


b(t)~B(f\ 
u«)~U(f ), 


(5) 


w{t)+->W{f). 

Applying the convolution theorem (Bracewell, 1978), 

B(f)=U(f)*W(fl ( 6 ) 

where the asterisk denotes convolution. In performing a 
straightforward autocorrelation estimation, one chooses be- 
tween the biased estimate in equation (1) and the unbiased 
estimate of equation (2). The advantages of the former are that 
the estimate satisfies the nonnegative definite property and 
that it deemphasizes the presumably less reliable values at 
larger lags. Its disadvantage is the same deemphasis, i.e., the 
bias in the estimate. The unbiased estimate eliminates this bias 
but sacrifices the nonnegative definite property. The method 
proposed in this paper eliminates the need to choose between 
these estimates by removing the bias, subject to satisfaction of 
the nonnegative definite property; it simultaneously offers esti- 
mation capability beyond either of the above approaches. 


The problem is to get an estimate of U(f) when B(f ) and 
W(f) are known. W(f) is the function 

W{f) — 1/N [sin 2 7rN//sin 2 */], 

which for large N becomes 

W{J ) ^ sin 2 nNf/Nn 2 f 2 = N sine 2 Nf 

Therefore, from equation (6), B(f) is a smoothed version of 
U(f). One way to remove this smoothing effect is to decon- 
volve B{ j). However, deconvolution is not straightforward if 
we wish to include the nonnegative definite constraint im- 
posed on autocorrelations, which corresponds to a nonnegati- 
vity constraint on U(f) (Papoulis, 1962; Robinson, 1980). Also 
note that deconvolving in the frequency domain produces u(t), 
which could be computed directly if the nonnegative definite 
property is ignored. 

In this paper a method to mitigate the effect of W(f) on the 
spectral estimate is proposed through a constrained iterative 
deconvolution process carried out in the spectral domain. 


ITERATIVE DECONVOLUTION 

The method of deconvolution using successive substitution 
(Bracewell and Roberts, 1954; loup, 1968; Lacoste, 1982; loup 
and loup, 1983; Jansson, 1984) was originally described in the 
time domain . For a function /(f), input to a linear shift in- 
variant system with impulse response h(t), the output x(f) is 

* 0 ) = J"/ (r')h(l ~ t ') dt' =f(t)* h(t), ( 7 ) 

where the asterisk denotes the convolution operation. Alter- 
natively, 


X(f) - F(f)H(n ( 8 ) 

where X(f\ F(f), and H(f ) arc the Fourier transforms of x(f), 
fUh and h(t) y respectively. The problem of deconvolution is to 
find /(/), given .x(t) and h(t). The principal solution (Bracewell 
and Roberts, 1954) for F(f) is 


f(/l , {«/»/«(/> 


for/: H(f)* 0 

for/: //(/) = 0. 


(9) 


loup and loup (1983) discuss the case when H(f) = 0. The van 
Cittert solution to the problem is stated in the time domain by 
the following equations: 


foU) = *(<)> 

/,(') =/o(') + 1*(0 -/o(t) * Mol. 

J (10) 


/.w =/,-.(o + -/„-!«) * m , 

where f n (t) is an approximation to f{t) which converges to the 
principal solution in a finite number of iterations. Time- 
domain constraints can be enforced after each iteration. loup 
and loup (1983) and the references cited therein discuss this 
method in greater detail, and Jansson (1984) gives a very com- 
plete development. 
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APPLICATION TO AUTOCORRELATION 
ESTIMATION 


The van Cittert solution is give ) in the time domain, 
whereas spectral estimation require . deconvolution in the 
frequency domain. In equation (6), B{ f ) and W(/) are known, 
and an estimate of U(f) must be computed that is positive at 
all frequencies. From equation (3), 

u{t) = b(t)/w(ty (11) 

The iterative method yields a solution lor U(f) given by 


u,(f) = u 0 (f) + |Vn - u, i/)* mf) , 


( 12 ) 


V m {f)=U ._ ,(/) + 


B(J) 




To understand these equations and t< * determine whether they 
satisfy conditions for convergence, we examine the time- 
domain equivalents: 


w o (0 = MO 

«i(0 = MO + 


b(t) - u { (f)w(f) 




(13) 


for which vv(f) is small or zero. These are the values of the 
unbiased estimate w(r) which are least reliable. 

Also, since u(0) and h( 0) represent the area of the power 
spectral estimates of V(f) and B(f) and since w(0) = 1, equa- 
tion (15} implies that the method preserves area. It is not 
preserved exactly due to the effects of constraints, truncation, 
etc. 


PRACTICAL CONSIDERATIONS 

An important consideration in iterative deconvolution is the 
growth of the length of the solution with each iteration, since 
the length of u n is longer than u n _ x by m — 1 when there are 
m samples in w(f). In practice it is necessary to truncate or 
limit the maximum length of u. For many response functions, 
only one expansion by convolution is necessary, and the effect 
of ignoring length expansion in further iterations is negligible 
( loup 1968; Hill, 1973). In some cases it may be necessary to 
use the neighboring replications of the power spectrum and 
the window' transform natural to the discrete Fourier trans- 
form representation. A program to implement this scheme was 
written and sample functions were tested. 

The pow r er spectral estimate is sometimes calculated directly 
from the square modulus of the discrete Fourier transform 
(DFT) of the data. Our method may still be applied to decon- 
volve the power spectrum, but only if the data are extended by 
at least (N — 1) zeros before calculating the DPT to avoid 
wraparound in the autocorrelation calculation. 


u n (t) = u„M0 + 

Upon substituting successively, one obtains 

«„w=MO{i+[i-*v(o]+r.i-H<o] J • •+[i-w(0]"}. (Hi 

This simple geometric series, which n ay be summed to give 

«„(» = { 1 — 1 1 - wU)]" M }MO/w<0, (15) 

converges for [ 1 — vv(f)| < 1. The detinition of w(f) guarantees 
convergence for all real input sequences except at / > /V, 
where w(N) - 0; and the discussion for zeros in the transfer 
function applies (loup and loup. 1983 }. In the limit as n tends 
to infinity, the convergent geometric cries sums to 

u(f) = b(t)/w{,). (16) 

This is precisely the result of the deconvolution [equation 
(11)]. However, note that for finite n the deconvolution is not 
perfect but is modified by the van < ittert window [equation 

(IS)]. 

This general form of iterative deconvolution is what is used 
for improved autocorrelation and p iwer spectral estimation. 
One key additional aspect, however, is the application of con- 
straints to equation (12). This is accomplished in an ad hoc 
fashion in that the nonnegativity is simply enforced on each 
U n {f) by setting all negative values to zero before beginning 
the (n + l)th iteration. The nonnegativity of the spectral esti- 
mate is thus ensured. 

In addition to the guaranteed no ^negativity, the iterations 
proceed gradually and can be terminated before too much 
bias is removed if the large-lag autocorrelation values are un- 
reliable. The constraints mainly affect those values of u„(0 at t 


EXAMPLES 

The following examples arc included to demonstrate the 
technique. The input sequence was extended by 32 zeros to 
allov* for autocorrelation extension. A 64 point FLT was used 
in all examples. Half the points in the autocorrelation and 
FFT represent negative times or frequencies and are not dis- 
played. The sampling rate may be considered to be unity, 
making the Nyquist frequency 0.5 Hz, and the plots are nor- 
malized so that the peak amplitude shown is the maximum 
ordinate value. The actual peak amplitude observed in the 
data is indicated on the ordinate of each plot. The increase in 
peak height which accompanies increased spectral resolution 
is therefore not readily apparent on the displays although it is 
indicated numerically on the vertical axes. 

Ten iterations were done in each example using the pro- 
posed algorithm. The number of iterations needed was not 
optimized in any way. The noise sequence added to the signal 
was approximately white. It was generated by transforming to 
time a spectrum of constant amplitude with random phase. 
Since only a part of this sequence was used, it is not strictly 
white. The signal-to-noise ratio was set to approximately unity 
in the examples with noise. 

The maximum entropy spectrum was obtained using the 
method described by Press et al. (1986). The order of the 
maximum entropy estimate was set to 5, with no attempt 
made to optimize the order. 

The first example is the 20 sample input of Figure l ; 

Input = sin (87it/32) 4- noise. 

The signal-to-noise ratio (S/N) was approximately unity. The 
steep triangular weighting on the autocorrelation due to the 
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Fig. 1. Single sinusoid with noise, (a) Original data, (b) Conventional autocorrelation, (c) Power spectrum [FT of (b)l 
(d) Autocorrelation estimate [inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolution from (c) 
(0 Maximum entropy power spectrum of (a). 


decreased number of terms at higher lag values is apparent. 
The autocorrelation obtained by iterative deconvolution 
decays more slowly. The extended part of the autocorrelation 
seems reasonable in that there is no dramatic change in its 
characteristics at long lags. 

Since the first sample in the spectrum represents /= 0, the 
spectrum should peak at the ninth sample, as occurs in Figure 
lc, the conventional spectrum. The iterative method places the 
peak at the eighth sample but the ninth sample amplitude is 
very close to the peak amplitude. The MEM spectrum places 
the peak at the correct sample. The peak amplitude in the 
MEM spectrum is smaller than that obtained by the iterative 
method. Side-lobe suppression is better in the MEM spec- 
trum. 

As a severe test of the proposed method, a nine-sample 
sequence was input in Figure 2. The input is sin (8 tu/ 32), a 
single period of a sinusoid with no noise. Due to the frequency 
and sampling rate, of the nine sample values, three are exactly 
zero. The autocorrelation is extended in a reasonable manner 
to four times its original length, and the spectrum is sharpened 
by the iterative method. The method gives reasonable results 
to lags of more than one and one-half times the original auto- 
correlation length. At subsequent lags, autocorrelation values 
are small. 

Again, one would expect the spectrum to peak at the ninth 
sample. However both the conventional spectrum and the iter- 


ative method peak at the eighth sample. The MEM estimate 
peaks at the seventh sample. This example was very sensitive 
to the order of the MEM estimate. Increasing the order from 5 
to 7 yielded totally erroneous results. 

In Figure 3 the 20 samples comprising the input are ob- 
tained from 

Input = sin (0.7 4- lnt/32) -f sin (0.9 + 9nt/32). 

In Figures 4a-4f, the same data and calculations are shown 
for the noise-added case. S/N was once again set to unity. Due 
to the presence of noise, the expected beating is not very evi- 
dent in the noisy input Figure 4a, compared to the noise-free 
input of Figure 3a. Comparison of the noisy autocorrelation 
Figure 4b to the noise-free autocorrelation Figure 3b, for 20 
points of this input, shows that the noise produces mainly 
modest changes in the autocorrelation amplitudes, except that 
the first negative lobe is less than half as large as it should be. 
In addition, noise causes the fourth (at a half-period), fifth, and 
sixth zero crossings to be delayed. The autocorrelations ob- 
tained by iterative deconvolution, Figures 3d and 4d, on the 
other hand, resemble much more the autocorrelation of 34 
samples of the original function, shown in Figure 3h, with a 
reduction in amplitudes at larger lags. The noise-free iterative 
deconvolution result in Figure 3e has all zero crossings in 
agreement with Figure 3h out to the 28th lag, except for the 
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Fig. 2. Single sinusoidal wavelet, (a) Original data, (b) Conventional autocorrelation, (c) Power spectrum [FT of (b)]. 
(d) Autocorrelation estimate [inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolution from (c). 
(f) Maximum entropy power spectrum of (a). 
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Fig. 3. Two sinusoids, (a) Original data, (b) Conventional autocorrelation, (c) Power spectrum [FT of (b)]. (d) 
Autocorrelation estimate [inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolution from (c). (f) Maximum 
entropy power spectrum of (a), (g) Autocorrelation calculated from (f). (h) First 32 autocorrelation lags using 34 
samples of the two sinusoids. 










Autocorrelation Estimation by CISD 


387 



FREQUENCY FREQUENCY 


Fig. 4. Two sinusoids with noise, (a) Original data, (b) Conventional autocorrelation, (c) Power spectrum [FT of (b)]. 
(d) Autocorrelation estimate inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolution from (c). 
(f) Maximum entropy power spectrum of (a). 
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Fig. 5. Square wave example, (a) Original data, (b) Conventional autocorrelation, (c) Power spectrum [FT of (b)]. (d) 
Autocorrelation estimate [inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolution from (c). (f) Maximum 
entropy power spectrum of (a), (g) Autocorrelation calculated from (f). 
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Fig. 6. Deconvolved prestack seismic data example, (a) Original data (b) Conventional autocorrelation, (c) Power 
spectrum [FT of (b)]. (d) Autocorrelation estimate [inverse FT of (e)]. (e) Spectrum calculated by iterative deconvolu- 
tion from (c). (f) Maximum entropy power spectrum of (a), (g) Autocorrelation calculated from (f). 
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sixth and seventh zero crossings which are slightly delayed. In 
Figure 4e, the noise has caused the following changes from the 
noise-free Figure 3e result. There are modest changes in rela- 
tive amplitude at several points and major changes in the first 
negative lobe, which is approximately one-third as large as it 
should be relative to the peak, and in the third negative lobe 
which is more than twice as large as it should be. The zero 
crossings agree, except for the third, which is advanced. In the 
noise-free case one may say that the autocorrelation of 20 lags 
of the function has been extended to good approximation to 
the first 28 lags of the autocorrelation of 34 lags of the func- 
tion by the deconvolution. Even when large amounts of noise 
are present, the extention can still be useful. The maximum 
entropy noise-free result, Figure 3g, does not resemble any of 
the extended autocorrelations considered by the authors, up 
to 40 sample points of input data. 

Using the conventional spectrum the two sinusoids are 
better resolved in the noisy data (Figure 4c) than in the noise- 
free data (Figure 3c), but this is an accidental effect of the 
noise. The noise-free iterative deconvolution spectrum (Figure 
3c) partially resolves two peaks of approximately the same 
height, Furthef resolution could have been obtained by using 
more iterations. The maximum entropy spectrum in Figure 3f 
resolves the peaks even better, but the peak heights and areas 
are very different even though the input sinusoids are of equal 
strength. For the noisy spectra, the iterative technique resolves 
the peaks completely (Figure 4e) with the same number of 
iterations as in Figure 3e. Again this is accidental and due to 
noise. The important point is that the noise in the convention- 
al spectrum (Figure 4c) has not been amplified relative to the 
main peaks by the deconvolution process. In the maximum 
entropy spectrum for the noisy data, Figure 4f, there is only 
one peak and the two sinusoids have not been resolved. 

A 20 sample symmetric square wave with a 4 time sample 
period is considered in Figure 5. The conventional auto- 
correlation is triangular. The autocorrelation obtained by iter- 
ative deconvolution has been extended in a manner that more 
closely resembles the autocorrelation of an infinite length 
square wave; there is some distortion apparent at large lags. 
In order to preserve the characteristics of the triangular wave 
autocorrelation, frequencies other than the dominant or fun- 
damental frequency, represented by sample 17 on the plots, 
must be present in the spectrum. This is seen in both the 
conventional and the modified spectra. The MEM method has 
resolved the dominant frequency, but other frequencies and 
harmonics are suppressed. The resulting autocorrelation is not 
quite as triangular as one would expect. In this example, 
MEM is perhaps not the most appropriate method. 

An example of real seismic data is presented in Figure 6. 
The input consists of the 20 samples of deconvolved prestack 
data; the peak of the conventional spectrum occurs at the fifth 
sample. It is apparent that although very few low-frequency 
components are present in the conventional spectrum, the 
maximum entropy method estimates a sizeable low-frequency 
component. Our method does not have this disadvantage. 

Since bias removal should increase autocorrelation ampli- 
tudes, we expect increased amplitudes for the iterative decon- 
volution (Figure 6e) and maximum entropy (Figure 6g) 
estimates compared to the conventional result (Figure 6b). 
Although both estimates show increases, the amplitudes of 
Figure 6e are greater, suggesting the iterative deconvolution 


approach performed better. Also, both methods extend the 
autocorrelations in a reasonably consistent fashion. 


EVALUATION 

Sometimes the iterative deconvolution spectrum is not as 
peaked as the one obtained from maximum entropy. The con- 
straints that lead to autocorrelation extension have the most 
effect at spectral valleys or notches that approach the zero 
baseline of the spectrum. If the spectrum that is input to the 
proposed method is flat and does not have notches in the 
spectrum, the constraints would be less effective. 

However, the experience of the authors indicates that the 
iterative method of autocorrelation and spectral estimation 
has the following advantages over the maximum entropy 
method: The iterative method does not assume an autoregres- 
sive model for the input data as does the MEM. The iterative 
technique proposed is fairly insensitive to noise. This insen- 
sitivity results from the technique’s attempting to improve the 
conventionally obtained spectrum through a gradual and 
stable iterative process. The maximum entropy method can 
result in a very different spectrum upon the addition of noise 
to the input. The maximum entropy method is sensitive to the 
order of the estimate. To do a reasonable job, some criteria 
(like Akaike’s criteria) must be applied. This implies running 
diagnostics prior to computing the final MEM estimate. The 
iterative method could be applied in a hands-off fashion to 
sharpen the conventional spectral estimate. The peak fre- 
quency in the maximum entropy method can err by one 
sample and may have incorrect relative and absolute ampli- 
tudes. The iterative method appears to suffer from this prob- 
lem to a lesser extent. 


CONCLUSIONS 

The iterative method cannot extend the autocorrelation in- 
definitely. The method’s ability to extend the autocorrelation 
results from nonnegativity constraints on the spectrum. Since 
the unconstrained iterative deconvolution is a linear process, 
it cannot add new lag values to the autocorrelation. 

In estimating the autocorrelation with the iterative method, 
the effect of the window used to produce the unbiased esti- 
mate is reduced, and the window is gradually made rectangu- 
lar as the iterations proceed. The autocorrelation is extended 
as necessary in order to maintain a nonnegative spectrum at 
the end of every iteration. The smaller values of the auto- 
correlation corresponding to the smaller values of the triangu- 
lar window can also be affected by the nonnegativity con- 
straint more than the large values. 

Naturally, there are no right or wrong answers when esti- 
mating power spectra (or autocorrelations) from short data 
sequences. Rather, the recorded data are assumed to be a 
segment of a larger sequence, and the ultimate issue is whether 
the model which generates the estimate is the same as that 
which generated the data. 

Finally, since the method is recommended primarily for 
short data sequences, computer cost is not a major consider- 
ation. 
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ABSTRACT 

Information contained in data that are in the form of a series of more-or-less resolved 
peaks is often unobtainable due to the limitations in resolution or response of the instru- 
ment. Adjusting the instrumental operating parameters to increase resolution usually has 
the effect of also decreasing the sensitivity and the signal-to-noise ratio, making detec- 
tion of small signals difficult. If a mathematical representation of a shift invariant 
instrument response function that describes the broadening effect on the input can be 
defined, then it is possible by deconvolution to restore the resolution to some degree.This 
process is represented by the solution of the convolution integral, which is achieved for 
many common data types through the use of iterative or Fourier transformation tech- 
niques. Although deconvolution techniques are becoming widely used, particularly in 
spectroscopic, acoustic, astronomical, and geophysical applications, and appear to be 
straightforward, care must be exercised to prevent the generation of spurious peaks 
which may be interpreted as being real data. This is particularly true when the higher 
frequencies in the Fourier transform are important in recovering the information. This 
paper describes both iterative and Fourier techniques developed in the course of on-going 
studies of the deconvolution process, and discusses some of the pitfalls which should be 
avoided. Results of recent research in optimizing iterative techniques and instrumenta- 
tion for deconvolution applications, and for evaluating and optimizing the efficacy of 
different methods of deconvolution for detecting peaks in specific classes of noisy data 
are also discussed. 

INTRODUCTION 

When instruments are used in the analytical sense, the data are typically obtained by 
periodically sampling the magnitude of the dependent variable. The independent varia- 
ble is most often time, but may be frequency, position, mass, wavelength, wave number. 




or any other parameter as well. The sensitivity of such an instrument is defined to be 
the smallest increase in the intensity of the signal representing the magnitude of the 
dependent variable that can be measured, with the detection limit being the smallest 
measurable signal that appears above instrument background and noise. The resolution 
or resolving power is the smallest increment in the independent variable that can be 
identified in the output and is therefore a measure of the ability to detect changes in the 
parameter being measured. 

The resolving power for a particular instrument is determined by the frequency response 
of the detection and recording circuitry and other factors influencing the output signal 
that represents a change in input. This is referred to as instrumental broadening, but in 
addition to the instrument response, it actually includes any parameters that affect the 
overall response of the system as well. Examples of these other sources are the inter- 
face between a sensor and the system upon which the measurement is being performed, 
mechanical or electronic limitations on the rate at which the independent variable can be 
scanned, and the presence of fundamental or environmental noise. Fundamental noise is 
that arising from the physics of the measurement process itself, while environmental 
noise results from externally imposed influences. With care, the latter can be made 
arbitrarily small, but fundamental noise has a limit below which it cannot be further 
reduced, therefore imposing a lower limit on detectability. The situation is, therefore, 
that a change in the input parameter is represented by an output signal whose response is 
dictated by the effects of instrumental broadening. Whenever these changes are slow or 
when only first order approximations are required, the signal obtained will represent 
adequately the variations in input. However, in many cases, valuable information is not 
observed or obtainable due to the broadening effects. 


Signal averaging or differentiation are both useful techniques to improve the resolving 
power of the instrument. The most effective computational method is deconvolution 
which, however, is mathematically difficult to apply in that it tends to amplify noise 
and has other difficulties. In some cases deconvolution will result in spurious peaks that 
may be interpreted as being real if not carefully examined. When properly applied, 
however, deconvolution of the signal can yield an improved resolution, often much great- 
er than can be obtained by careful tuning of the instrument or otherwise enhancing the 
signal, and sometimes even exceeding the theoretical limit of the resolving power. The 
improvement in resolution is determined in part by how accurately the response function 
(impulse response) can be determined for each point at which the deconvolution is to be 
carried out, and how completely the effect of noise can be addressed. 


If the instrument response does not change appreciably while the measurement is being 
carried out, the response is called "shift, invariant," and the relationship between the 
observed signal h(x) and the input (or ideal) signal f(y) is defined by the convolution 
integral 

+ ® 


h(x) 


f(y) g(x-y) dy - f * g 


( 1 ) 
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where g is the response function representing the broadening effects. The discrete 
form of this equation is given by 
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Complete discussions of deconvolution and reviews of many of the techniques are 
given by Frieden,* Robinson and Treitel, 2 < Mendel, 29 and Jansson. 22 

In this paper a background in basic Fourier deconvolution is summarized, and mention 
is made of newer Fourier transform techniques which improve upon it. The source of 
sidelobes is discussed as well as why deconvolution generally amplifies noise and is 
therefore ill-conditioned. A review of iterative deconvolution techniques is given 
including recent developments. Examples of the pitfalls due to sidelobes and noise are 
presented. The optimization of iterative techniques is discussed as well as a new 

method for evaluating different deconvolution methods for their ability to detect true 
peaks and reject false peaks. 

FOURIER DECONVOLUTION 

Fourier deconvolution proceeds directly from the Convolution Theorem.* This theo- 
rem holds that if h(x) is the convolution of f(x) and g(x) as in Equations (I) or (2), 

then H(s) the Fourier transform of h(x), is the product of the Fourier transforms of 
g(x) and f(x): 


H(s) - F(s) G(s) 


( 3 ) 


As long as G(s) is unequal to zero the transform of the ideal function f(x) can be 
recovered by division, but when the experiment fails to transmit frequencies that are 
present in the ideal function, that is, G(s) is zero for frequencies s for which F(s) is 
not zero, information about f(x) contained in these frequencies is lost, and f(x) can not 
be perfectly restored. Bracewell and Roberts* suggest that an initial approach in such 
cases might be to define a principal solution F (s), an approximation to the transform 
of the ideal function, whose value is zero whenever G(s) is zero. Thus, 


F p (s) - 


H(s)/G(s) 

0 


(s.C(s) * 0) 
{ s : G ( s ) - 0} 


( 4 ) 


The approximation to the ideal function, f (x), is then obtained by taking the inverse 
ranstorm of F p (s) t and the process by which this resolution enhanced solution is 
obtained is called simple inverse filtering. 


The loss of information about F(s) at frequencies for which G(s)=0 often creates 
unwanted characteristics in the approximation f (x). In many experiments the trans- 
it"? the T esponse f unctI on is non-zero (except perhaps at a few isolated points) for 
all frequencies below a certain critical "cut-off" frequency s c and zero above this 
frequency. In these cases F p (s) is truncated. This usually gives rise to spurious peaks 
(Gibbs oscillations) which are not present in f(x) and greatly complicates the interpreta- 
tion of spectra enhanced by simple inverse filtering. 

As the Convolution Theorem clearly shows, the transform of the measured function 
H(s), should be zero at values of s for which G(s) is zero. Noise at these frequencies 
cannot have been transmitted through the system and is therefore called incompatible 




noise. 30 Noise at frequencies s for which G(s) is not zero is called compatible noise. 
Incompatible noise can be removed by bandpass filtering without further loss of infor- 
mation before Fourier deconvolution if G(s) is known. Compatible noise which is 
present at frequencies for which G(s) is small is an additional serious obstacle to accu- 
rate resolution enhancement by simple inverse filtering. Filtering which reduces or 
removes compatible noise necessarily forfeits part of the signal, and information 
needed for the accurate restoration of f(x) is lost. In particular, noise tends to be more 
important at high frequencies so that low-pass filtering which decreases the bandpass 
of the approximation to the transform of the ideal function broadens the function it- 
self, thereby partially tending to defeat the purpose of resolution enhancement. 

For additive noise, the principal solution has the form F (s) = [H(s) + N(s)]/G(s) 
where N(s) is the transform of the noise. The example shown in Figure 1 includes 
plots of sample |H(s)+N(s)| and |G(s)| functions and the resulting |F_|, which is the noisy 
F(s) truncated at the cutoff of G(s). As the cutoff frequency is approached from 
below and G(s) becomes smaller, the magnitude of N may exceed that of H, and hence 
noise predominates and is amplified by the deconvolution. In this case, the high 
frequency component of the transform domain spectrum is often deleted, a process 
defined as "low-pass inverse filtering." This reduces the effect of the noise but does 
not eliminate the Gibbs oscillations and it decreases the resolution as discussed above. 
(The application of a taper (window) to the transform can reduce these oscillations. 38 ) 



j iSUre |L\ ° f lh * Fourier tranjfom ’ ° { th « principal solution for the transform of noisy measured 

data |H(s) + N(s)| is the magnitude of the transform of the noisy data. The majnitude of the transform of the 
impulse reponae is labeled |G(s)|. Simple inverse filtering yields a transform with magnitude |F (s)|. 

A number of techniques have been developed to improve inverse filtering by restoring 
the lost resolution and correcting its tendency to create spurious peaks while controlling 
the effects of noise. Many of these methods use constraints to incorporate special 
knowledge of the ideal function into the deconvolution process. Deconvolution con- 
strained to produce an approximation to the ideal function f(x) that has as little nega- 
tivity as possible, 12 13 14 frequently referred to as either Minimum Negativity or the 
Howard extrapolation, is an example of this approach to deconvolution. A similar 
method was also developed by Gerchberg. 9 Many experiments produce measurements 
of intensities which can be either positive or zero, but for which negative values can 
be created only by unwanted noise or sidelobes from deconvolution. Howard uses the 
minimum negativity constraint to extend the frequency spectrum of a deconvolved 
function expected to be positive by retaining initially only a small number of the 
lowest Fourier coefficients and restoring the lost coefficients one at a time so that each 




successive coefficient produces an approximation to f(x) that minimizes the sum of the 
squares of the negative values of the previous approximation. 

The minimum negativity constraint has been applied to microwave spectra by Howard, 
extending the spectrum of h(x) that has been curtailed by the measuring spectrometer. 
L. Wang and Rayborn 38 have extended the principle of minimum negativity directly to 
the deconvolution process itself, forming F p (s) and using the principle of minimum 
negativity to extend F p (s) beyond the cut-off frequency by first making f (x) an even 
function. p 

Another technique for restoring the spectrum of a function truncated in deconvolution 
is to approximate the high frequency spectrum which has been truncated in deconvolu- 
tion with the high frequency spectrum of a function of the expected shape. 41 11 This 
method for extending the bandwidth uses simple inverse digital filtering to establish 
the size and location of peaks in the output of a laboratory spectrometer. The ideal 
peak shape for the instrument is determined either fiom a theoretical understanding of 
its performance or by experimental study and observation of an isolated peak. An 
artificial function is then formed by superposition of peaks of ideal shape of the sizes 
and at the locations determined by the inverse filtering. The Fourier transform of the 
artificial function is taken and the high frequency portion of this spectrum is extracted 
and used to replace the high frequency portion of the spectrum of the inverse filtered 
function which was truncated either naturally by the low bandpass of the spectrometer 
or deliberately because of the presence of noise at the high frequencies. Extending the 
Fourier spectrum in deconvolution by splicing the high frequency spectrum of an arti- 
ficial function has been found to decrease the size of sidelobes created by inverse 
digital filtering by up to about 50 percent. 

ITERATIVE DECONVOLUTION 

One of the commonly used approaches to deconvolution is an iterative one, first pro- 
posed in its simplest form by van Cittert in 1931. 35 Excellent reviews are given by 
Frieden 8 and Jansson. 22 G. Ioup and J. Ioup 17 summarize additional literature. The 
form of the van Cittert iterations is as follows: 

- h + (h - h*g) 
f 2 “ fi + (b - f^*g) 

(5) 


f n " f n-l + < h - f n-l*S> • 

The gradual nature of the iterations causes a controlled simultaneous increase in resolu- 
tion and noise, allowing a compromise to be determined by the selection of iteration 
number. Because the iterations are accomplished in the function or time domain (as 
opposed to the Fourier transform domain), function domain constraints can be applied 
easily in an ad hoc fashion at each iteration. If the signal-to-noise ratio of the exper- 
imental data is low and more noise control is needed than that provided by the van 
Cittert iterations, an iterative noise removal technique proposed by 
Morrison, 30 15 17 2S 26 applied before the deconvolution, can give additional noise 
control. 




The convergence conditions for the Morrison and van Cittert iterations are well estab- 
lished. 8 1S 10 11 * 17 There are many analytically and experimentally determined 
response functions for which the iterations do not converge. A modification of the 
deconvolution iteration called the reblurring or mirror image procedure which con- 
verges for any initial response function has been developed by Kawata and Ichioka 23 
and independently by LaCoste. 24 Jansson 22 gives background information and a 
discussion. An alternative always-convergent iterative noise removal and deconvolution 
technique was given by G. Ioup 18 and applied to two-dimensional data by Whitehorn 39 
and Whitehorn and G. Ioup. 40 

One of the principal objections to the use of the iterative approach to deconvolution is 
’the fact that it can be very slow for long data sets and impulse responses or for wide 
impulse responses. To overcome this problem the research group at the University of 
New Orleans has been investigating accelerated iterations and single filter application 
of the iterations in the Fourier transform domain. 3 31 The single filter application is 
based on the transform domain representation of the iterations. The last of Equations 
(5) may be written in the transform domain as: 

F n - F n-1 + <» - F n-1 G > • < 6 > 

By successive substitutions one may solve for F n in terms of G and H to obtain 

F n - [1 - (1 G) n+1 ]H/G . (7) 

Similar results have been obtained for Morrison’s noise removal 15 and the newer 
convergent iterative techniques of noise removal and deconvolution. 39 40 18 Use of the 
so-called van Cittert or equivalent window makes possible the accomplishment of many 
iterations as a single filter. Because multiple convolutions imply expansions of the 
duration of the solution function in its independent variable, there is a possibility for 
serious wraparound error 32 Amini et al. 3 and Ni et al. 31 have shown that for many 
experimental data types wraparound error is negligible. 

PITFALLS 

Examples of effects on deconvolution of Gibbs oscillations and noise clarify the diffi- 
culties of the deconvolution process. Figure 2 contains the results of two approaches to 
deconvolution for mass spectrometric analysis of a gas containing oxygen and 
methane. 21 Simple inverse filtering (not shown) gives large sidelobes which are the 
Gibbs oscillations. The irregular nature of these oscillations is due to the interference 
of the sidelobe patterns of the two main peaks. The artificial function approach 
(Curve B) reduces the spurious peaks somewhat but does not eliminate them. Iterative 
noise removal followed by iterative deconvolution with a nonnegativity constraint 
included gives Curve A. Because deconvolution with a normalized impulse response 
should preserve areas, the elimination of the negative spurious lobes by the use of a 
constraint reduces the positive lobes as well. 21 Since the interference of the Gibbs 
oscillations due to the presence of two main peaks causes the adjacent positive and 
negative sidelobe areas to differ from each other, the area cancellation is incomplete 
and two small positive lobes remain after iterative deconvolution at mass to charge 
16.0559 and 16.0963. These might be interpreted as small additional mass peaks if one 
is not aware of this pitfall and the need for careful analysis. 
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Figures 3 and 4 show the effects of noise on a deconvolution of the synthetic data 
given by the curve in Figure 5, which is generated by convolving an ideal function 
consisting of three narrow-width, separated Gaussians of heights 100, 100, and 50 with 
a medium width Gaussian impulse response. For each example (Figures 3 and 4) 
different Gaussian distributed function domain noise has been added at a signal-to- 
noise ratio (signal peak to noise standard deviation) of 1000 to 1. Even though such 
noise is too small to be seen on a graph, the noise amplification by deconvolution (ill- 
conditioning) gives large noise spikes in the result of Figures 3 and 4. The original 
input function is also shown for comparison. If a wide Gaussian impulse response had 
been used the noise iinpiificntion would have been even worse 


One striking feature of Figures 3 and 4 is how much the deconvolved results differ 
from each other. Differences in the noise become dramatic after deconvolution. This 
behavior casts doubt on the ability of deconvolution methods to be optimized based on 
an average noise spectrum. Certainly any such optimization should be tested on specif- 
ic noise realizations to understand whether the use of average noise properties is 
appropriate. An alternative approach, discussed in this paper, is to do many noise 
cases in a simulation and then to calculate statistics. 
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Figure 3. Deconvolution of the synthetic broadened data h of Figure 5. Unmarked curve is original input. 
Curve marked with x'j is deconvolution of a noisy h with SNR=1000. Noise if Gaussian distributed. 
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Figure A. Same as Figure 3 except that a different noise set (still Gaussian distributed) is used to produce a 
noisy h with SNR=1000. 



Figure 5. Input f and output h after convolution of f with a medium width Gaussian. 




OPTIMIZATION OF ITERATIVE DECONVOLUTION 

A long unanswered question for iterative methods has been how many iterations to use 
for a given data set. Experience and some trial and error for each data type have 
provided crude answers, but a systematic investigation did not begin until ten years 
ago. Wright 42 initially investigated requirements for optimization of Morrison’s itera- 
tive noise removal. 30 18 Her work was extended by Wright and G. Ioup 43 and J. Ioup 
and G. Ioup. 10 Finally a very complete analysis was performed by Leclere 25 and co- 
workers. 20 Since then the method has been extended to van Cittert’s iterations and to 
Morrison’s and van Cittert’s iterations applied sequentially. 21 The first optimization of 
noise removal and deconvolution iterations used sequentially was for G. Ioup’s always- 
convergent method by Amini 1 and co-workers. 2 4 The reblurring procedure has also 
been optimized. 2 

The initial hypothesis for the optimization studies is that for a given data type and 
signal-to-noise ratio (SNR), simulation with many different noise sets, and testing each 
for optimum iteration numbers, leads to reasonable statistics for the iteration number. 
It is assumed further that by varying the SNR, a smoothly varying set of the mean 
iteration numbers versus SNR can be generated with standard deviations, as well as 
that these standard deviations will not be too large. Although many problems present 
themselves in the course of the work, the final result is that useful curves of average 
iteration number versus SNR and average mean squared error of the deconvolved re- 
sults versus SNR can be generated. Studies include narrow to wide Gaussian impulse 
responses as well as synthetic exploration seismic data. Several different noise types 
are included over SNR’s from less than 1 to over 1000. 

To answer the question of how changing the input model would affect the results, an 
input consisting of a rectangle function followed by a triangle function is substituted 
by El Saba 7 for the three- Gaussian input shown in Figure 5. The optimization gives 
different optimum mean iteration number versus SNR curves, as expected. It is there- 
fore important for experimentalists to do an optimization for the data type of interest. 
While there is some computer time involved, the labor is minimal since the optimiza- 
tion code exists. Once optimization results are obtained for a general data type, no 
further simulations are needed for experiments of the same type. 

A by-product of the optimization of noise removal and deconvolution iterations used 
sequentially is the first solid answer to the question of when iterative noise removal 
was helpful and when it was not needed. By plotting the mean squared error after 
deconvolution versus SNR for sequential-use optimization on the same graph as the 
result for the optimization of deconvolution iterations alone, one can decide below 
which SNR noise removal is needed. The mean squared error will be smaller with 
noise removal at the lowest SNR’s, but above some SNR value the noise removal will 
not improve the result significantly or at all. The mean squared error curves for a 
fairly narrow Gaussian impulse response are shown in Figure 6 and for a wide Gaus- 
sian in Figure 7 using the input given in Figure 5. The narrow Gaussian deconvolution 
is improved by noise removal up to a SNR of approximately 90. The wide Gaussian 
deconvolution is improved by noise removal for all SNR’s included in Figure 7, i.e., up 
to SNR 150. 

Although the mean squared error (L2 norm) is principally used in this work, the LI 
norm for optimization has also been tested. It is important to emphasize that t is 
method is not limited in the choice of optimization criterion, and that many others 
could be used. We are currently studying non-norm type measures. 
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Figure 6. Mean squared error versus natural logarithm of SNR (a) after optimization of iterative deconvolution 
alone (solid line) and (b) alter optimiration of sequential application of iterative noise removal and deconvolu- 
tion (dashed line) for a fairly narrow Gaussian impulse response. 
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Figure 7. Same as Figure 6 with (a) iterative deconvolution optimzied alone (dashed line) and (b) iterative noise 
removal and iterative deconvolution optimized together (solid line) for a wide Guassian impulse response. 

These studies have also provided a methodology for deciding how to optimize instru- 
ment design and operating parameters to achieve the best deconvolution results. 19 In 
experiments with a compromise between SNR and resolution, an approach should be 
available to decide hew to configure the instrument to obtain the needed experimental 
data. If deconvolution of the data is to be part of the process, then the optimization 
approach should include the deconvolution. Anuni et al. 4 have shown that a three- 
dimensional plot can be created with a surface of mean squared error after optimiza- 
tion of iterative deconvolution plotted versus SNR and resolution (impulse response 
width, for example) as the independent variables for systems with one degree of free- 
dom. See Figure 8. Once the curve of SNR versus resolution is established for a given 
instrument in the SNR-resoiution plane, an upward projection of this curve will give 
another curve at the intersection with the surface. 19 The instrument resolution value 
corresponding to the minimum in the latter curve should give the best deconvolved 
result. 
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Figure 8. Surface of mean squared error after optimiration of iterative deconvolution versus SNR and Gaussian 
inpuise response width. The latter is a measure of resolution. 

EVALUATION OF DECONVOLUTION ALGORITHMS 

Systematic evaluation of deconvolution algorithms seldom takes place. Particular 
algorithms sometimes are applied to output of particular instruments and the results 
compared with known results of the instrument at high resolution in a qualitative 
fashion (e.g., J. Ioup et al. 51 ). The technique may be used or discarded depending on 
whether it has distorted the spectrum unacceptably, created false, or ghost peaks, or 
simply used too much computer time. Sometimes synthetic spectra are generated, 
convolved with an assumed instrument function, and used to test the deconvolution 
algorithm by determining the square of the difference between the original spectrum 
and the deconvolution As a figure of merit for evaluating deconvolution algorithms, 
this square difference criterion produces a large value for a spectrum which is faithful- 
ly deconvolved if it is merely translated slightly. Recently, the problem of evaluation 
of resolution enhancement has been attacked by regarding it as a problem of detecting 
an unresolved peak in the presence of both noise and the larger, fused peak. 38 37 33 
This point of view seems directly applicable to those spectral problems in which the 
presence or absence of a peak is the primary question. It also permits techniques long 
applied to the detection problem to be applied without modification to the problem of 
resolution enhancement, permitting for the first time judgements as to which deconvo- 
lution techniques are superior when the principal objective is to find a weak, fused 
peak and which are superior when the objective is to avoid falsely believing a weak 
peak to be present. Receiver Operating Characteristic, or ROC, curves are plots of the 
probability of falsely detecting an artifact of the deconvolution. Their application to 
the resolution problem shows that the efficacy of ideal low-pass filtering before 
deconvolution of a spectrum of fused Lorentzian peaks depends, of course, on the 
bandpass of the filter. For the example shown in Figure 9 a superior deconvolution as 
judged by the ROC curve produced occurs when about eleven of the 128 Fourier 
coefficients are retained. Of particular interest, however, is the fact that filtering with 
a narrower or wider bandpass, although each produce inferior deconvolutions, affects 
detection in different ways. Narrow-band filtering produces deconvolutions which are 
superior in avoiding false detections at low probabilities of detecting true peaks, while 
wide band filtering produces deconvolutions which better detect signal peaks when a 
relatively high rate of false alarms can be tolerated. When the gross efficacy of low- 
pass filtering prior to the application of an inverse digital filtering deconvolution algo- 
rithm is evaluated by maximizing the area under the resulting ROC curve, optimum 




deconvolution is found to occur when only ten Fourier coefficients are retained in the 
low-pass filtering. When the efficacy of the filtering and deconvolution is evaluated 
by minimizing the square of the difference between the deconvolution and the original 
spectrum,, a gentle minimum and, hence, optimum deconvolution is found when nine- 
teen coefficients are kept. Thus, optimum deconvolution in the sense of minimum 
square difference does not produce optimum ability to detect peaks without creating 
spurious ghost peaks. Furthermore, when deconvolution with a minimum negativity 
constraint was evaluated using ROC curves 38 it was found that, although minimum 
negativity produced enhanced spectra whose peaks were easily interpreted visually, the 
ability to detect small, fused peaks in noise by means of a matched correlation filter 
was only marginally superior to optimum low-pass filtering. This result casts doubt on 
the deconvolution algorithms. It is clear that much work remains to be done in im- 
proving evaluation of deconvolution techniques; nonetheless, viewing resolution en- 
hancement as a problem in detection offers promises of quantitative evaluation of 
deconvolution with a numerical appreciation of the risk of misidentifying a small peak 
discovered in deconvolution of noisy, fused spectra. Numerical assessment of the 
presence of ghosts and the failure to detect true peaks should greatly increase the utili- 
ty of all deconvolution algorithms. 
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Figure 9. ROC curves (probability of detection versus probability of false alarm) for low pass filtering by 
varying the threshold at signal-to noise ratio 3.183. 'Square: cut-off number 39; triangle: 12; diamond: 3; cross: 
no cut off. 
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(2) how much analog precision is needed in the connections in the net- 
work, (3) the number of training examples the network must see before 
it can be expected to form reliable generalizations, and (4) the efficiency 
with which a network extracts information from the training data- 

John Denker, Daniel Schwartz, Ben Wittner, Sara Solla, John flop- 
field, Richard Howard, and Lawrence Jackel, Complex Systems, in press 
(1987). 


L Q1 8 2J An Analog VLSI System for Neural Network 

^ESgfflmg Experiments DANIEL B. SCHWARTZ and RICHARD E. 
HOWARD AT&T Bell Laboratories . 

Because the complexity available from standard VLSI has grown far be- 
yond our ability to simulate it, it has become interesting in its own right. 
Adaptive neural network models are an example of a class of complex 
systems where a mapping directly onto VLSI is of great practical and 
fundamental interest. However, the continuously variable connections 
required for adaption are not easily represented in a digital world. We 
are building a collection of analog circuits from standard digital CMOS 
with variable strength analog connections ba*ed upon charge storage by 
a pair of MOS capacitors. The capacitors are tied together by a string 
of FETs, allowing the connection strength to be monotonically varied by 
moving charge between them. Our current designs have 7 bits of ana- 
log depth of both polarities. The chips have about 10 2 3 connections and 
can easily be cascaded to make larger networks. The available computar 
tional speed is dominated by i/o bandwidth of the host controller. We 
will discuss use of such chips and their limitations. 
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018 3 Jgpwjfic Beal for * Boren, gjatemwiiii Anh« - 
^PfUiwslxY . J4.S. Wartak, C.Y. Fong, Department of 
Physics, University of California, Davis. — We used the 
model Hamiltonian 
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LOU 5 J Surface Loss In s Parabolic-Equation Model . 
'wAWW* E.M. HEAD and V. JOBST, Naval Oceanographic 
Office and ELEANOR &, HOLMES, Science Applications 
International Corporation. — Ocean-surface lose of 
acouetlc energy la often given as a function SL(e) of 
the grazing angle e. If p(z) la the coaplex acoustic- 
pressure field (froa a parabolic-equation aodel) as a 
function of depth z near the surface, a Fourier trans- 
fora F(I) of p(z) yields pressure as a function of the 
vertical wave number L E la proportional to sin #, 
thus F(C) la a function G(e) of •• We account for the 
surface loss by multiplying G(e) by a lose function U» 
— related to SL(e) — before transforming back to 
physical apace. The method also is applicable to bottom 
lose. Numerical Implementation, angular resolution, and 
limitations of the method are discussed. Numerical 
examples are presented. 


#30 1 

^ $18 6 ^Affects of Noise fin Pressure and Modal Amplitude 
Matched Field Processors . * GEORGE M. FRICHTER, IV, 
JULIETTE V. IOUP, Unlv. fif fifiit Orleans .** GEORGE B. 
SMITH, Xavier Unlv. . GEORGE E. IOUP,** Pnlv. fif New 
Orleans . CHRISTOPHER FEUILLADE , Svntek . GRAYSON H. 
RAYBORN, Unlv. al Southern Miss. .** snd DONALD DEL 
BALZO, NORDA - -Modal amplitude matched field processing 
for acoustic signals received by a vertical array of 
hydrophones Is used to determine the effects of 
spatially correlated and uncorrelated noise fields on 
pressure and modal amplitude matched field processors. 
Various amounts of white isotropic noise snd spatially 
correlated noise as calculated by a normal mode noise 
model are combined with the field due to e submerged 
acoustic source to produce simulated cross spectral 
matrices A phona- to-aode spacs mapping is then used to 
obtain the corresponding cross amplitude correlation 
matrices Both conventional and maximum likelihood 
processing are used. Results show that spatially 
uncorrelated noise degrades modal amplitude processors 
more then spatially correlated noise. 

★★Supported In part by ONR/NORDA Contract N00014-87 -K-600 


to study the thermodynamic properties of the one- 
dimensional boson system with on-site anharmonicity, 
and with A much smaller than c. For the calculation 
of partition function we have used the path-integral 
method. The Dyson equation is solved in the nearest- 
neighbor approximation. The resulting expression for 


the free energy is evaluated in the static approximation 
using the steepest descent method. The behavior of spe- 
cific heat for different values of f and A is examined. 




Color Induced Transitions in the Presence of a Nonlinear 
Potential. GTpT TSIR0N1S. P. 0R1G0L1N1. University of 
California, San Diego.- We show that the negative diffusion 
coefficients exhibited by current approaches to the Fokker- 
Planck equation for non-Markovian and bistable processes 
result from .assuming that the system reaches a conventional 
steady state 1 . By lifting this assumption* 4 we show that when 
a critical value of the noise correlation lime t is exceeded, 
the process of escape over a barrier agrees with an exact 
prediction for the large-t regime and thus that the linear 
response approximation behind our theory produces exact 
results for arbitrary correlation times. 

, » 

1. G. P. Tsironis, P. Grigolini, "Rate processes activated by 
color noise: Bridging two exact limits’, UCSD preprint 

2. J. Masoliver, 0*. J. West, K. Lindenberg, Phys. Rev. A 

35, 3086 (1987) 


JOSEPH E. MURPHY. University of New Orleans New Orictni. 
LA-70148; STANLEY a7c^-BINO ^ 
and Development Activity. NSTL. MS 39329-5004. - 
Underwater acoustics is usually not discussed at APS meetings, 
but rather is confined to peer review meetings. However, given 
the close proximity of the Navy's lead ocean environmental 
RDT&E laboratory, the Naval Ocean Research and Development 
Activity (NORDA) located 45 miles from New Orleans, we take 
this opportunity to present a review of ocean acoustic propagation 
modeling. In the ocean, the index of refraction is variable; 
acoustic transmission paths are curved and the coupling of the 
refracted* reflected, and diffracted acoustic fields from 
boundaries give rise to complicated classical physics problems. 
The prominent acoustic mod els are based -oo perm^l pa ode, 
parabolic approximation, FFT, and modified ray methods., JJach 
of these include a limited number of physical mechanisms. We 
have therefore developed a coupled full-wave range-dependent 
ocean acoustic propagation and scattering model based on the 
finite element method. This model is superior especially at low 
frequencies. Numerical simulations will be presented showing the 
effect of a fractal under-ice interface with ice keel on the fully 
coupled range-dependent underwater acoustic field. 
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Orleans **- -An important design and parameter selection 
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, h od which combine* a modification of earlier suggestion* 1 
■> j . u b«pace diagonal! ration, realising advantage* of both meth- 
5- for a matrix expressed in a planewave ba*i*. Comparison 
f MT m P totic convergence rate, for several iteration method. 
t hat the combined method i« a', least a factor of three 
’itter than the best previously known method. The method 
^implemented in a planewave basis with separable nonlocal 
oseudopotentials so that computational effort scale* a* NlogN 
* bwi* of N planewaves. 

. wr Payne, J.D. Joannopoulos, D.( Allan, M.P. Teter, and 
1 pH. Vanderbilt, Phy,. Rev. Lett. 56, 2656 (1986). 
j A Williams and J. Soler, Bull. Am. Phys. Soc 32 , 562 
(1987). 


A Uji 2 , S olution to the Initial Value Problem for the 
Quantum Nonlinear Schrodinger Equati on. M. J. POTASEK and 
mM /: bYURKE, AT&T Bell Laboratories - -The quantum nonlinear 
WU Schrodinger equation provides an integrable quantum field theory 
Spit- that has been solved by a number of methods. Most recently, 
Gutkin 1 - 2 has developed an intertwining operator technique for 
£Kdr‘- obtaining the time evolution of the field operators. Using 
Gutkin’s formalism, we show how to obtain the exact time 
?*( evolution of an initial Schrodinger state vector. The suitability 
't- of this formalism for numerical computation with application to 
pulse propagation in optical fibers will be discussed. 

fPt;. j. E. Gutkin, Ann. 1st Henri Poincare 285 (1986). 

2. E. Gutkin. J. Func. Anal. (1987). 


glv.n l»puls. re.pons* and data type, above which It 1. 
not beneficial to precede the deconvolution with the 
nolle removal . 

♦♦Supported in part by NASA Grant. NAG-1-485 and HAO- 1-804 
1(5 g P ioup. Bull. An. Phys. Soc. 26, 1213 (1981) 


w iv 185 i single Filter Application af Always -Cgnvgr & 3Ht- 
Ite fatiye neconvolutlon ifl Physical Patft a haihong ni, 
ABOLFAZL M. AMINI, TAHAR a BENSUEID, GEORGE E. IOUP, and 
JULIETTE W. IOUP, Univ. lifiw -Single filter 

application of any iterative technique, when possible, 
presents significant computational economic advantage, 
but it should not be used until Its performance is 
evaluated against that of the iterations. Wraparound 
errors associated with a finite length DPT calculation 
of the filter must be considered. The optimization of 
the always -convergent iterative technique of Ioup for 
noisy data is reported In the preceeding abstract. In 
this investigation the sensitivity to wraparound of the 
DFT single filter equivalent window is established by 
gradually increasing the zero padding of the data for 
peak to standard deviation signal - to - noise ratios 
varying from 10 to 150. It is found that the wraparound 
error is small enough to be negligible, even when almost 
no zero padding is used. These results show that very 
rapid application of Iterative deconvolution to physical 
data is possible. 

♦★Supported in part by NASA Grants NAC-1-485 and NAG-1-804 
*G, E. Ioup, Bull. Am. Fhys. Soc. 26, 1213 (1981) 


N18^y Statistical Optimization ol Morrlaa ii-I Itcrat l- Y tt 
Removal and van Clttcrt'a UtlflUYA EtCPimlut l PIU 
JAKES H. LECLERE, GEORGE E. IOUT, and JULIETTE V. IOUP, 
Pnlv. fil Orlsani **- -Korrison* s iterative noise 

removal and van Cittert's Iterative deconvolution have 
been used to remove experimental broadening for various 
physical data typaa. 1 Heretofore the number of 
iterations needed and other conditions of use have been 
determined by trial and error. We have developed a 
atatistical optimization procedure to determine optimum 
use of the methods for any computer- generated noise type 
and any signal -to -noise ratio domain of interest. We 
report the results for LI and L2 norm optimization and 
several noise types for a signal -to-noi.e ratio domain 
from 2 to ovar 1000. The contrast between point 
successive and point simultaneous iterations is also 
discussed as is the effect of allowing the deconvolved 
result to expand as the iterations proceed. 
Combined optimization of the two techniques is 
presented. 

♦★Supported in part by NASA Grant NAG-1-485 
V E. Ioup and J. W. Ioup, Geophysics 48, 1287-1290 

(1983) 


N 18 4/ Qpti nt zation Cpnvargcnt I terat i n g 

^ ^opvolutiOD fax Physical ABOLFAZL M. AHINI , 

GEORGE E. IOUP, and JULIETTE W. IOUP, VnlY -*. HfiJ* 

Orleans **- -Statistical computer simulation is used to 
optimize the always -convergent iterative noise removal 
and deconvolution technique of Ioup. 1 By considering 
the mean square error versus iteration number for 50 
noisy data sets', one can calculate the mean optimum 
iteration number and its standard deviation, a a well as 
the average mean square error and its standard 
deviation. D®ta with peak- to -standard- deviation signal- 
to-noiae. ratios (SNR) varying from 10 to 150 are 
considered. By applying the iterative deconvolution 
alone, it shown that there exists an SNR, for any 

Vol. 33, No. 3 (1988) \ 
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N186^j Hon-S.pnr.bl. Huia.rlc. 1 8v«lu|>tlon of 
VGV-Wa tn fch. Schwlnq.r Variational wthod . 
C.A. Weathertord, Florida A6M U.j* A numeri- 
cal method is presented for the evaluation of 
the VGV term which appears in the denominator 
of the Schwinger variational expression for 
the T-matrix.** 2 The method employs an evalua- 
tion of a partial differential equation for the 
GV part, and then is followed by the calcula- 
tion of a two dimensional integral. An applica- 
tion to electron scattering from a minimal 
basis set H2 model is presented. The possibil- 
ity for efficient evaluation on vector compu- 

‘Supporte^by^SF grant PHY-8711805 and NASA 
grant NCC 2-492 

1. W.M. Huo, T.L. Gibson, H.A.P. Lima, and V. 
McKoy, Phys. Rev. A36, 1632 (1987). 

2. W.M. Huo, M.A.P. Lima, T.L. Gibson, and V. 
McKoy, Phys. Rev. A36 , 1642 (1987) . 

12:12 — - 

N18 7 Generalized^ F ourier series f or non-li near quantum 
mechanics.- J. DIAZ BEJARAN0 and A. MARTIN SANCHEZ, 
Unlversid ad de Extremadura, Badajoz . — A simple 
generalization of the usual Fourier series uBing the 
generalized exponential and circular functions is 
presented. The functions themselves are developed in 
a new, more simple way. They are solutions of common 
linear and non-linear wave equations. The series are 
given in terms of Jacobi elliptic functions in a form 
as similar as possible to the usual Fourier 
presentation. Several examples are given that correspond 
to the most usual textbook Fourier series. 

Thenks are due to CA1CYT (project n» 1179-84). 
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An Ex tension Theory for Lattice Green F unctions. 
j[ $CHWALrt and HVk/SCMWALK, U. of N. Dakota--A method 
©resented of using the known Green functions or den- 
ties of states (DOS) for a given lattice Hamiltonian H 
YJnd 6reen functions or DOS for any lattice Hamllto- 
fan 11® expressible as a rational function of H, l.e. H e 
xOO- The ^ orma H sm I s further developed to permit u- 
i-jg the known Green functions and DOS of two Hamllton- 
** H and K to obtain those of any Hamiltonian H e in the 
ra generated by the direct products (Hal) and (I«K), 
I Is the Identity. Namely, 

l a^,, (H m iK n ), with 3,^, real numbers. 

^r?lts may be obtained either analytically or numerical - 
< Other properties such as electrical conductivity are 
%o extended from H and K to H ? . A study Is presented 
two examples, (1) The Cartesian product of two Sler- 
1©sk1 fractal lattice Hamiltonians, also a hierarchical 
>ictal which Is infinitely connected with spectral di- 
mension 4 ln3 / in5 * 2.730- • • (11 ) 2D and 3D plaid lattice 
^{fcmlltonlans formed as products of ID Fibonacci chains. 
^Jhese lattices are quasl-perlodlc and admit Inflation/ 
^deflation transformations, but do not have 5-fold rotation. 

CL.; ) 


Si 


If 10 Rydberg States of the Rare-Gas Van der Waals 
\mr% . WlN6 Y1 PU and t. H. GREENrr Lou1s1ana~Sute 
C*LHilt1 channel quantum defect theory Is adapted to 
describe bound and autolonlzlng Rydberg states of the 
fire gas dimers. As In our earlier paper* , related to a 
similar study by de Prunele*, a nonperturbatlve Ferml- 
type analysis combines readily with MQOT, giving an 
imtormous simplification. This permits the description 
complicated avoided crossings among Rydberg state 
ptentlal curves* e.g. for ArHe, XeNe, without requiring 
large-scale ab Initio calculation. Autolonlzlng 
Resonance structures In the photoionization cross 
lection are also calculated between the fine-structure- 
split ionization thresholds, accounting partially for ^ 
observations of Dehmer and Pratt* 9 
^Supported in part by the National Science Foundation 
"1. Y. Du and C. H. Greene, Phys. Rev. A 36, 971 (1987) • 
\ de Prunele, Phys. Rev. A 35, 496 (1987); also Phys. 

A 36, 3015 (1967). 

M. Behmer and S. t. Pratt, J. Che*. Phys. 77 , 4804 
(1962). ~ 


SESSION 119: GENERAL MECHANICAL 
PROPERTIES AND NDE 
^Wednesday morning, 23 March 1986 
Mardl Gras K at 8:00; 

IL C Cook, presiding 

SOO 

119 1 Acoustic A xma In TI. HI AndRI. RII Laua Group 

grystali, DAVID T. CHUNG, Howard University. Thera 

are tlx and seven elastic constants for TI (Rl) and Til 
^ (RII) group .crystal* respectively. In an earlier paper, 
Khatkevich 1 indicated that the crystals of TI (RI) are 
related to TII (RII) by a rotation of an angle y about 
the 4- fold (3- fold) axis. This rotation of acoustic 
axis is the only distinction between the two groups so 
far as the acoustic properties are concerned. In the 
present work, we like to show that by the use of so- 
called invariant constants, this rotation ot <f> comes 
out naturally from the inherit properties of these cons- 
tants. Invariant constants 2 * are the elastic constants 
which are independent of the specific coordinate system 
belr^t used. The detail expressions of <f> tor TI .TII and 
RI, RII groups will be presented at the meeting. 

1. A. G. Khatkevich, Sov. Phys. Crystallography, 6,361 
(1962). 

2. T. P. Srinfbasan, J. of Math, and Mech. , 19, 1019 

(1970). 


8:12 

119 2 The Frenkel-Kontorova Model With Nonconvex 
Interparticle Interactions and Strain Gradient^ S Mananer 
and A. R. Bishop. LANL .-- We study the statics and dynamics of 
a chain of atoms moving in a periodic potential with nonlinear, 
nonconvex internarticle interactions, and with strain gradients 
which we model by including next nearest neighbors* 
interactions through the discrete Hamiltonian H = 
L n u n 2 /2 +a(u n + i-u n j 4 -B(u ri+ i-ii n / 2 +y(u n + i-2u n + u n .i) 2 — cosu n . 
We obtain the phase diagram within an ansatz of periodically 
modulated configurations. These generalize the homogeneous 
(for P < 1/8) and dimerized (for p>i/8) configurations already 
reported for Y = 0, and are given by: u n = na\ + 6, for n-l...N, 
and u ft = fia 2 + f >2 for n = jV + 1...N + M. The dynamics of 
tranistions between different configurations when the 
parameters are varied is also investigated and we show that 
these are dominated by nucleation processes, which occur on 
short time scales compared with the subsequent slow growth. 
Possible relation of the model to the dynamics of twin 
boundaries recently observed in the copper-oxide high- 
temperature superconductors is discussed. 

(***\ 

\I19 3/ Tempe rature Dependence ..of Third Order Elastic 
amg(tAXita _Q f. ,KMpF 3a * w. CAO, G.R. BARSCH, Penn State 
LU, W. JIANG, M. A . BREAZEALE, 11 + of Tenneaaee . --We have 
measured the three nonlinearity parameters along the 
principal synvnetry directions for KMnF^ from 298 to 
350K by means of acoustic second harmonic generation. 

In conjunction with our earlier data on the temperature 
dependence of the pressure derivatives of the elastic 
constants the complete set of the six third order 
elastic constnts has been determined in this 
temperature range. For c lllf c 112' c 123 and c 166 
temperature dependence is linear, indicating that the 
effect of the ferroelastic transition at 186K (manifest 
in elastic anomalies) Is no longer present above 300K, 
and permitting us to eliminate the effect of zero point 
and thermal motion by extrapolation to absolute zero. 

The static T.O.E. constants thus obtained differ 
significantly from the R.T. values. Both static and 
R.T. values exhibit large deviations from the Cauchy 
relations. The results are also compared with those 
for other perovskites . 

•Supported by Office of Naval Research under Contract 
No. NC0014-82-K-0339. 

4^6 

119 4 ANflVa-Convtmnt Iterative Deconvolution f£i 
AtOUttU Non -Detractive Evaluation. EDVARD J. MURPHY, 
JULIETTE V. IOUP, and CEORGE E. IOUP, Unfv af BfiSt 
Orleans . DORON KISHONI , Coll . of William and Mary and 
MSA Langley Res . Cen . - -Acoustic energy sources 
generally have a finite time duration and a ringing 
shape which can make the evaluation of individual 
reflections difficult. Deconvolution can be an 
important tool for signal analysis. In this work the 
detected acoustic signals are deconvolved using the 
Always - Convergent Method of Ioup. The Always - 
Convergent technique Is applied to data recorded during 
the quantitative analysis of materials through Non- 
Destructive Evaluation In which ultrasonic signals are 
used to detect flaws in substances such as composites. 
An Important part of processing the signal la the 
normal J zation since it is useful to know the strengths 
of the reflections. Various methods of normalization 
are investigated and the most effective method is found 
to be the one which uses the change in the sum of the 
absolute values of the amplitudes in the signal before 
and after processing. Results of the application to 
actual data are shown. 

1 G. E. Ioup, Bull. Am. Phys. Soc . 26, 1213 (1981) 


r8:48 

119 Thermal Properties of PEEK (poly-ether- 
etTier-Vetone) Based Materials * Ryan Giedd, Ewan 
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-howcd that ihe directivity index (L)1 j ol jn intimieiy densely populated 
lineal shell array was about equal to that of a sphere. Extrapolating to 
^screte elements, this means the shell requires far fewer elements. They 
did not compute DI using amplitude shading, due to the impractical cost 
0 f such systems at that time. Today's technology removes that constraint, 
'phis work revisits the problem with shading, using an approach for choos- 
. the amplitude shading coefficents tha= maximizes D! [H.S.C. Wang, 
j Acoust.Soc. Am. 57, 1076-1084 ( 1975} j. Calculations have been made 
for the DI of shaded cubic volumetric arrays, forming beams perpendicu- 
lar to one of its faces, in the presence of is< iropic noise. Results show that 
for 27 and 125 element arrays with element matrix spacings of 1/2 wave- 
length. a full 10 log ( number of elements) tan be obtained for DI Work is 
underway to investigate larger arrays and >malier spacing*,. The approach 
w j|| also be extended to nonisotropic no se fields. [Work supported by 
flORDA and NOSC exploratory develoj ment programs ] 


edge oi individual element amplitude and phase response tor accurate 

resu Its. Two in situ methods of array calibration are described and results 
from the September experiment are presented. The first method used 
transmissions from a low-frequency source of known location and power 
level. Simulating the conditions encountered during the transmission, the 
power arnving at the array was predicted by several acoustic propagation 
models. By comparing the array response at specific frequencies to the 
response predicted by the models, an absolute calibration was obtained. 
An error curve for the phase data was generated by unwrapping the phase, 
accounting for a sampling offset in the array, and subtracting a multiple 
linear regression curve. The second method determines relative amplitude 
levels by examining the average ambient noise power output of a specified 
frequency band across the array. Using spectral, coherence, and direction- 
alit) plots, the level of self-noise in the array was shown to be below that of 
the ambient noise being measured. These two independent methods pro- 
vide a consistent set of element calibration values used for array beam- 
forming. [Work supported by ONT ] 
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9:00 

H7. Least-squares and single-filter always-convergent iterative 
deconvolution of transient signals for correlation processing. James 
H Leclere, George E. Ioup, - ’ Juliette W. loup, -> and Robert L. Field 
(Code 244, NORDA, Stennis Space Center, MS 39529) 

Correlation processing for distributed sensors is most accurate for 
short pulses and those whose autocorrelation is sharply spiked. For longer 
transient signals, multipath arrivals at ea :h sensor have significant inter- 
ference with each other, and it is difficult to identify individual arrival 
times. Deconvolution of the received signal to sharpen the transients is 
one method to decrease the overlap and increase the accuracy with which 
travel times can be identified. Deconvolution can also be applied after 
cross correlation to sharpen the autocon elation of the transients. Least- 
squares deconvolution is the most commonly used approach for acoustic 
signals. It has the disadvantage of being computer intensive when filters 
for long transients are needed. An alternative approach, the single-filter 
application of the always-convergent iterative technique, is faster and pro- 
vides variable control for noise The two techniques are compared for 
actual underwater acoustic multipath transient signals. Single filter appli- 
cation of always-convergent iterative noise removal is compared to the use 
of a modified Blackman-Harris window for noise control. -l Also at the 
Department of Physics, University of New Orleans. 


9:05 

H8. Comparison of double and triple cross correlation for arrival time 
Identification of amplitude- and frequency-modulated acoustic transient 
signals. Juliette W. Ioup, -1 George E Ioup, -1 Robert L. Field, and 
James H. Leclere (Code 244, NORDA, Stennis Space Center, MS 
39529) 

The triple cross correlation of three signals is a simultaneous function 
of two lags. It is an alternative to cross correlations taken two at a time for 
determining the lags for a given sourer at three distributed sensors. It 
should offer improvement in arrival time identification only when the 
statistics of the signal have significant third moment components. In this 
study, amplitude- and frequency-modulated snythetic transient signals 
are propagated over several possible paths to three sensors, and the triple 
correlation of the received pulses computed, as well as the cross correla- 
tions of the same three signals two at n time. The efficacy of these two 
approaches is compared for a variety of amplitude- and frequency-modu- 
lated transient signals and multipath interference conditions. Also at 
the Department of Physics, University < if New Orleans. 


9:10 

H9. In situ acoustic calibration for « large aperture array. Barbara 
J. Sotirin (Marine physical Laboratory A-005, Scripps Institution of 
Oceanography, La Jolla, CA 92093) 

During September 1987, a large ape t ure acoustic array was deployed 
vertically in the Northeast Pacific to study low-frequency noise in the 

S17 J. Acoust. Soc. Am. Suppi 1 , Vol. 84, Fall 1 988 


9:16 

H10. Abstract withdrawn. 


9:20 

H11. Matched-mode processing corrections for array tilt and bottom 
type. James A. Mercer (Applied Physics Laboratory, University of 
Washington, Seattle, WA 98105) 

In a related effort, Homer Bucker has shown that matched-mode pro- 
cessing for an unknown sound-speed environment can be significantly 
improved if correction factors for the mode-line amplitude functions can 
be determined. The correction factors are obtained when a source with 
known location is available to calibrate the system. This paper describes 
the results of applying the same techniques for simulated cases of un- 
known array till and bottom characteristics. 


9:25 

Hi 2. Self-consistent modeling of signal and noise in a three-dimensional 
environment. John S. Perkins, W. A. Kuperman, and F Ingenito (U S. 
Naval Research Laboratory, Code 5160, Washington, DC 20375-5000) 

Previous propagation work is extended to model surface noise, ship- 
ping, and signal sources in a fully three-dimensional environment. The 
noise cross-spectral density matrix for a vertical array is computed as the 
sum of a local contribution and propagation from distant small patches of 
ocean surface. Propagation from any point to the array is made efficient 
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Role of Remanent Magnetization in the MAGSAT CtusiaI Anomaly 
Field in the SW Indian Ocean 

L Fulknon and J Roark (Astronomy Program, Univcrsuy of 
Maryland. College Park. MD 20742), 

H Frey and H Thomas (Laboratory for Terrestnxl Physics/Goddani 
Space Right Center. Greenbelt MD 20771; 301 286-5450} 

The MAGS AT regional crustal magnetic anomalies m the SW 
Indian Antarctic Ocean arc doe to a combination of induced plus 
viscous remanent magnetization and TRM in Cretaceous Quiet Zone 
(KQZ) regions. Two broad, roughly parallel, SW to NE trending, 
multi-peaked lobes of positive reduced -o- pole (RTP) anomalies 
dominate the region, one lying south of Africa and the other north of 
Antarctica Some of the peaks of these anomalies correlate well with 
the location of submarine platcuas which are tectonic conjugates, 
i.e., formed together but now separated. But the shape, location of 
many of the peaks and amplitude contrast of the northern lobe, 
which runs from the Agulhas Plateau northeastward to the 
Madagascar Ridge, appears to be controlled mostly by TRM in KQZ. 
crust structural characteristics (i.e., thickened crust) account for 
only about 20% of the tool anomaly amplitude. Based on modeling 
results, the TRM contribution varies from about 10 A/m over the 
Mozambique Plateau and Basin) to about 3 A/m over the Agulhas 
Plateau. Tmnskei Basin and Madagascar Ridge (TRM assumed 
distributed through layer 2). This inferred differential TRM is 
consistent with available drill core data. The southwestern portion of 
the Enderby Abyssal Plain has a 3nT positive MAGSAT anomaly 
over it which is centered south of the Conrad Rise. This entire area 
is KQZ crust, and the anomaly seems little related to the Conrad 
structure. The centroid of another 3 nT anomaly which lies between 
the Maud Rise and Astrid Ridge may also be controlled by the KQZ 
crust rather than the structural features which flank this portion of 
the extreme southwest Enderby Basin. Overall there is good 
correspondence between the MAGSAT RTP anomalies observed for 
conjugate plateaus and adjacent ocean basins in that portions formed 
together at the same time but now well separated seem to have 
similar TRM contributions to the local anomaly contrasts observed 
Surprisingly, it may be possible to use MAGSAT data to infer ihe 
limits of KQZ crust where this is poorly known. 


Magnetization Contrast of the Pacific Cretaceous Qulel Zone 
Based on Magsat Data 

P B Toft and I Arkam-Hamed (Dept, of Geological Sciences, McGill 
University. Montreal, Quebec. H3A 2A7, Canada; 514-398-8052) 

The absolute value of magnetization of oceanic lithosphere is poorK 
known. Information from drill holes is only for ihe uppermost 0.. 
km and that from ophioliles is scattered and sparse. Magnetn 
anomaly inversion gives only a value of magnetization contrast. 

A constraint on the absolute value may be obtained from Magsai dat. 
and the Pacific Cretaceous Quiet Zone (QZ) of normal polarity 
There are few Magsat anomalies over the east Pacific, and most larg- 
anomalies in the west Pacific are correlated with plateaus an.! 
seamounts. An anomaly over the Cretaceous Hess Rise, for example 
is modelled from topography and crustal thickness with a iota 
magnetization contrast of 10,000 A (magnetization i thickness). 

In the central Pacific, however, are anomalies that are not obvsousL 
correlated with topographic features. Some of these may be due I > 
a magnetization contrast between the QZ and its surroundings: if ih - 
magnetic signatures of the narrow bands of normal and revers-’ 
magnetizations surrounding the QZ sum to zero at satellite aliitudr. 
and if there is no susceptibility contrast across the QZ boundary, the i 
Magsat anomalies may result from an edge effect of the QZ. 

To lest this hypothesis, the magnetic anomaly of the QZ is calculate 1 
at 400 km and it is filtered to simulate removal of long wavelengths 
overlapping the core field that are extracted from Magsat data along 
with the core field. The QZ with a total magnetization of about 
10,000 A produces 500-1000 km wavelength features spatially assoc. 
ned with the QZ boundary, which are similar in magnitude, wav- 
length. and location to observed Magsai anomalies Both ihe Paaf i 
QZ and isolated plateaus, such as the Hess Rise, indicate a tot ! 
magnetization of about 10.000 A for the Pacific oceanic lithospher- 
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The Ifitemltr of Magnetization of Ocean k Basalts as a Function 
of Age; Compiled Data from ODP and DSDP Basalts. 


rock samples. There data were filtered to remove sample* with 
anomalous geochemistry and from non mid ocean ndge 
environ menu. There new results show that, on average. J for oceanic 
basalts (1) steadily decrease from 0 to 35 Ma, (2) increase from 35 to 
50 Ma, and (3) exhibit no significant change between X) and 140 Ma. 
The mean value of J all studied basalts is 3.5 A/m. 

The observed decrease in J from 0-35 Ma is consistent with 
our understanding of the effect of progressive k-w temperature 
oxidation on the intensity of natural remanent magnet ration. The 35- 
50 Ma increase mutt be due to either a submarine process which 
effects oceanic basal ts of this age. or a change in the intensity of the 
paleomagneoc field. The lack of variation in J betwren 50 and 140 
Mi is no / consistent with recent models which predict predici that 
basalts formed during the Cretaceous Normal Magnetic Superchron 
should have substantially higher values of J than basalts formed 
before, or lfter, the Superchron. Therefore, we suggt si that the high 
magnetic fields measured by MAGSAT over the Cretaceous Quiet 
Zones in the Atlantic ocean result from (l) a thicker extrusive layer, 
or (2) an increased contribution from the lower, intrusive layers of 
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Deconvolution for Increased Resolution in AIM 


Clyde J Bergeron, Jr, George 5 IdMg. 

Juliette W loup, Long B Trinh, and Abolfazl 
M Atnini (Physics Departaent and Geophysical 
Research Laboratory, University of Hew 
Orleans, New Orleans, LA 7014B; 

(504) 286-634 1) 

Deconvolution is a standard technique for 
removing the effect of the instrument or 
other response functions fro* data. In 
airborne electromagnetic (ACT) eeasurements, 
there is an effective impulse response for 
the A EH measuring device due to the large 
footprint of the device. We present 
approximate line and point impulse response 
functions calculated from the Modified I»ag® 
Method (HIM) representation of the A EH field. 
We apply these functions in an iterative 
deconvolution of data produced from two- 
dimensional sodels. The deconvolved results 
in qeneral show a increase in the effective 
resolution of the AEM data. (Work -upported 
in part by the U. S. Army Cold Regions 
Research and Engineering Laboratory.) 


i from ODP and DSDP Basalts. 
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RANDALL MACKIE and THEODORE MADDEN (Both «: 
Dept, of Earth, Atmospheric, and Planetary Scienc es, M I T., 
Cambridge. MA 02139) 

We have developed numerical algorithms for computing the 
electromagnetic response of a fully inhomogeneous, 3D earth model 
due to a uniform current source far above the earth (thu is *•** 
magnetotelluric response). Our algorithms arc finite difference 
algorithms, but they arc based on the integral forms of Maxwell s 
equations rather than the differential forms. This eliminates ®e need 
to approximate derivatives of earth properties; indeed o ne only needs 
be concerned with the issue of taking averages of earth properties 

Finite difference algorithms invariably lead to large systems of 
equations to be solved, especially for realistic 3D earth models. We 
have investigated relaxation methods (conjugate direction algorithms) 
and direct methods for solving these systems quickly and accurately 
The relaxation solutions are quick, give reasonable answers, and do 
not require large amounts of computer storage. We have found that a 
multiple scaling technique used in conjunction with relaxation 
methods works especially well. The direct solutions are more 
compuier intensive than the relaxation methods because they can 
require large amounts of storage space and involve doing many 
matnx inversions The solutions from our algorithms compare well 
with the solutions from Wannamaker's integral equation algorithm. 


J E Pima and H P Johnson (School of Oceanography, University of 
vtashmgtoo, Seattle. WA 98195; 206-543-8542) 

H. Sakai (DepL of Earth Sciences, Toyama University, Toyai na. 
Japan) 

In contrast to the original Vinc-Matthews model of he 
magnetization of oceanic crust, wc now know that oceanic rocks are 
subject to changing physical and chemical conditions which have j he 
potential to modify the magneffc properties of crustal rocks he 
variation in the intensity of magnetization (J) of oceanic basalts ^as 
initially investigated by Bleil and Petersen (1983), who comps if a 
paleomagnetic da ta from DSDP Legs 1-65. A more recent study at the 
University of Washington focussed on a re -ex ami nation of dri led 
basalts (0-155 Ma). and includes all of the DSDP and new ODP bird 
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HIM Inversion of AEM Data for Groundwater 

Juliette W loup, and Clyde J Bergeron, Jr 
(Physics Departaent and Geophysical 
Research Laboratory, University of Hew 
Orleans, New Orleanp, LA 70140; 

(504) 286-6715) 

The Modified laage Method (HIM) fer inversion 
of airborne electromagnetic (AEM) 
measurements can b« applied to groundwater 
studies. The effective depth and 
conductivity of the groundwater is determined 
from the simultaneous measurements of the AEM 
transmitter /receiver system altitude end the 
complex low frequency secondary field. The 
inversion of the high frequency A1TM data 
allows a determination of the average 


conductivity of the overburden layer. The 
results of auch an analysis of data from a 
groundwater survey performed by the Dighe* 
Company in Michigan will be presented. 
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Determination of Paleoenvironmental Conditions From 
Soils Using Rock Magnetic Techniques 

Michael J Singer (Dcpi of Land, Air and Water Resources, Univ. 
of Calif., Davis, CA 95616) 

Bruce Moskowiiz and Kenneth L Verosub (Both at: Dept, of 
Geology, Univ. of Calif., Davis. CA 95616) 

Pinchas Fine (Institute of Soils and Water, Volcani Center, P O 
Box 6, Bet Dagan. ISRAEL) 

Soils are sensitive indicators of paleoenvironmental conditions, and 
horizon by horizon chemical analysis of a soil can provide detailed 
information about climate and climate change. We have been using 
rock magnetic techniques to supplement traditional wet chemical 
methods for tracing the movement and transformation of iron in 
soils. Our work has shown that enhancement of the magnetic 
susceptibility of a roil is not a surficial process, as was previously 
believed, and that the increase in susceptibility results from the 
retention of inherited magnetite as well as from the precipitaion of 
pedogenic maghemite. Both of these factors are influenced by the 
water content of the soil, and the enhancement process ceases when 
a soil becomes poorly drained We have also found that 
morphological discontinuities found in soils often have a 
corresponding magnetic susceptibility anomaly so that susceptibility 
measuiments can be used to evaluate the suitability of soil sequences 
for detailed pedologicaJ analysis. Our one enigmatic result is that the 
ratio Xttarm •PP^rs to be independent of soil horizon and particle 
size fraction within a soil, This result implies that the ferromagnetic 
particles in a roil fall within a narrow size range 
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Diagnosis for Grslglt* (r*jS 4 ) io Cretacaoms a«d», 
Dorth Slops , Alaska » and Dolocsoo tedtMts, Lsk* 


R L Reynolds , A Nicholson, K I Coldhaber (USCS, HS 964, 

Sox 25046” Denver, CO 80223, 303-236-J 3u3 ) 

S N Colssn (USCS, Woods Hols, MA 02543) 

J U Ring ( On 1 v . of Rhode Island. Narragenaett , R1 

02882 ) 

C A Rice. M L Tuttle, end 0 H Shaman (USCS, Darwer, CO 

80225) 

The presence of pos tdepos It lonsl grelflce 
ferriaagnetlc) can distort daposltlonal .agnatic record* 
and thereby jeopardise paleoenvironmental 
Interpret at ions baaed on sagnetlc suecept 1 bl 11 ty (MS) 
and remanent sagnet i zee Ion. Out studies have provided 
aeans to Identify gretglte sod evaluate Its effects. 

Using X-ray diffraction, tnemooagnetlc signature, and 
Mbesbaoer spectra, we nave identified abundant grelglte 
in Upper Cretaceous Budatones fro« the Stepson 
Peninsula , North Slope, Alaakj , but have found little 
evidence of it in Holocene Bud frt» southern Lake 
Michigan. 

Urclglte doeinates magnetic properties of_the Cretaceous 
audstonea, which have a Bean MS of 5.9x10 (vol SI) and 
a Bean NKM magnitude of 6.6xlU" 2 A/b. Geochen lcally . 
these Budatones resemble many Recent aarlne sediments: 
Acid-volatile sulfur (AVS; from grelglte and nonmagnetic 
monosulf lde) ranges from 0.021 to 0.2UX (by weight), 
disulfide S (from pyrlte) ranges from 0.02X to 0.46T. 
and the ratio of AVS plus disulfide S to organic carbon 
averages 0 . 34 , 

comparison of magnetic properties with the distribution 
ol sulfur species in Holocene aud from southern Lake 
Michigan suggests the presence of grelglte and pyrit*. 
Neither mineral, however, appears to obscure the 
detrltal magnetic record, except perhaps In shallow {<20 
cm) intervals in some cores. In these shallow 
Intervals, three cores achieve maximum sulfur values 
(AVS U.02 to U.lbt; disulfide S U.G11 to 0.24 I.) mostly 
within the ranges of the Alaskan mudstones. In one 
core, the AVS profile mimics the MS profile, both 
reaching their maximum* (AVS 0.021; MS 4.4x10 ) at 6 

cm. In the upper lu cm ot another core, high AVS 
content (maxiaura 0.16Z) corresponds to relatively high 
MS (s.lxUT*) and NKM (l.lxlO -2 A/m). both of these 
results suggest that some of tne AVS may be in the form 
of grelglte. In the latter core, however, high MS also 
correlates with high sand content, indicating that 
detrlial oxides contribute muen more to MS than does 
grelglte. This finding is supported by examination ot 
subnet ic separal c 
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Applications of the Modified Image Method to the AEM Study of Sea Ice 

INTRODUCTION 

Classical image theory may provide a simple method for 

calculating the secondary field produced by a conducting 

structure in response to an applied field. The three-dimensional 

screening distribution that is induced in the conductor is 

replaced by an image (or a distribution of images) of the primary 

source. When the primary field is magnetic, classical image 

theory may be used in the case of a ''perfect” conductor (high 

freguencies and/or high conductivities) or for the response of a 

superconductor to a constant magnetic field. The general 

criterion for the applicability of classical image theory is that 

the weighted average depth of the screening current 

distribution (the superconducting penetration depth or the complex 

skin depth in the ohmic case) be much less than the altitude of 

the primary source above the conducting surface. 

The essential constraint on the image sources for the 

magnetic case is that they produce a secondary field at the 

surface of the conductor whose normal component, B , is equal 

sn 

but opposite to the normal component of the primary field. B 

' pn' 

thus satisfying the continuity condition on B n at the interface. 

A symmetry constraint on the tangential components, B . and B 

pt st' 

results in their equality. Thus 

B (surface) = *B^ + - 2 B t T , ( 1 ) 
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where T is a unit vector in the conducting surface. The 
conductor surface is generally termed the image surface. 

Only simple interface geometries allow for simple 
distributions of image sources. A plane interface requires only 
one image source to satisfy the electromagnetic boundary 
conditions on B n - In this case the secondary field produced by 
the induced current distribution is accurately replicated in the 
nonconducting half space by the field of a single image source. 

An alternate source distribution to the image source is provided 
by a surface current *K which is given by K = n x B at the image 
surface, where 3 is the normal to the surface. Thus K is given 

by 

t B T« (2) 

where T' is a unit vector tangent to the surface but 
perpendicular to T and n. The secondary field produced by the 
surface current distribution K is identical to the field produced 
by the image source, but this calculation is less direct since it 

involves an integration over the surface. 

The modified image method (MIM) has all the elements of 

the classical image theory EXCEPT that the image surface is 
relocated to one weighted average screening length BELOW the 
conducting surface. For the ohmic case this distance is complex 
and is given by [exp(-i*/4)/y2]5, where S is the electromagnetic 
skin depth 6 = J( 2/Ov**)- For a P lane layered conducting 
medium, 6 is modified by a complex correction factor, Q (Bergeron 
et al., 1987). The complex image field produced by this 
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assumption is in good agreement with the secondary field given by 
the one-dimensional Sommerfeld theory (Sommerfeld, 1909) . 

For non-planar conducting models (two- and three- 
dimensional) , there is no general prescription for determining a 
distribution of image sources that will satisfy the 
electromagnetic boundary condition, and in general none exist. 

The standard methods of calculating the secondary field produced 
by two- and three-dimensional structures (e.g., the finite 
element method (Lee and Morrison, 1985)) are generally of an 

iterative nature and hence computationally slow. 

HIM provides a fast, efficient, but approximate method 

for calculating the secondary field. We assume the general 
validity of Eg. (2), which defines a surface screening current 
distribution K in terms of the tangential component of B p . This 
screening distribution is at a complex skin depth, 
[exp(-W 4 )/ 72 ]«, below the conducting surface, where now n is^ 
the LOCAL normal to the conducting surface. Equation 2 gives K 
in terms of the primary field alone. The heart of this 
approximation is that the relation B gt - B pt is still valid. ^ 
In the context of the MIM theory, the primary field B p 
on the image surface is a formally complex function, since it is 
a function of the primary source strength (real), the lateral 
displacement of the source from the image surface point (real) , 
and the vertical separation between the source and the image 
surface point, h + 5 eff (complex) . It follows that and * in 

the image plane, given by Eg. (2), and the secondary field 
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generated at the detector by K are complex, i.e., the secondary 
field is not in phase with the primary field at the detector. 

Thus in the image-source surface a discontinuous change in a 
model parameter results in a local discontinuous vertical 
displacement of the source surface. That is, a discontinuous 
change in a y a ^ or d ^ for a two-layer model, results in a 
discontinuous change in S f , and hence a vertical shift in the 
source surface. 

In this approximation the source surface is now 
disjointed and the strength of the secondary surface distribution 
is calculated locally in various areal cells of the surface. The 
individual contributions of these source surface cells are summed 
up at the detector coil location. In general the area of a 
surface cell is decomposed into Cartesian components 

, A A A A 

da = da n = a da i + p da j + 7 da k , 

where a, p, and 7 are direction cosines with respect to the x, y, 
and z axes, respectively. 

There are two facts that render the assumption of the 
general validity of Eq. (2) at least plausible. The first is 
that the field is screened from the interior of the conductor 
independent of its geometry, i.e., B i nte rior “ °* The second 
is that Eq. (2) provides for the approximate satisfaction of the 
electromagnetic boundary condition at the image surface. *K 
produces a stepwise discontinuity in at the image surface from 
2 B pt to approximately zero, and likewise implies that B n = 0 at 
the surface, thus satisfying the boundary conditions on B n * 
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The secondary field produced by K of Eq. (2) is easily 
and rapidly calculated, taking about 1 sec/survey point for 
fairly arbitrary two-dimensional models. See Table 2. 

The ultimate utility of this approximate method of 
calculating secondary fields produced by two- and three- 
dimensional structures depends on how well these MIM fields agree 
with secondary fields generated by other methods of calculation 
and with scaled model laboratory measurements and survey results. 

In this report we give the resultant numerical 
calculations of the MIM secondary fields produced by models of 
ice keels. The specific models studied were: 

* 

1) rectangular trough models in which we examine the dependence 
of the secondary maximum (in ppm of the primary field) and the 
secondary field half maximum width (in meters) on keel depth, 
width, and bird altitude. 

2) a triangular ice keel model which we label the Berkeley 
model, identical to one used by Becker et al. (1987) . 

We compare the results of the MIM two-dimensional field with the 
Berkeley two-dimensional calculation and to a one-dimensional 
Sommerfeld calculation. 

3) a "CRREL" model which is based on an Arctic sea-ice survey 
ground truth data s;et provided by CRREL (bird altitude, ice and 
snow freeboard, and keel depth versus range) . 

THEORY 

In this Ejection we include the "working equations" 
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which are used in the numerical computations of the secondary 
fields produced by the current distribution given by Eg. (2) for 
two-dimensional structures. We also include the HIM and 
Sommerfeld one-dimensional formulations which result in fields 
that exactly track the geometry of the model. The two- 
dimensional MIM approximation produces fields that vary more 
smoothly but still generally track the model. There is one 
exception to this rule which will be pointed out and discussed 


later in this report* 

Table 1 illustrates and defines the notation used for 
the normalized secondary field detected by the various possible 
permutations of horizontal coil pairs. 


Table 1 


coil pair 

name 

horizontal coplanar 
horizontal coaxial 
horizontal mixed 1 
horizontal mixed 2 


diagram 

4 <=> 


T X 


G~ 0 

T X 


0 

T X 



normalized secondary field 

ZZ - H sz / H p 

XX » H SJ /H p 

ZX - H sx /H p 

XZ = H sz /H p 
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Table 2 gives the total CPU time per survey point for a line 
survey and the horizontal surface element dimensions. 

Table 2 


Field 


XX 

XZ 

ZZ 

ZX 


} 

} 


sec/survey point 

1.34 ax 


1.74 AX 


step sizes 


0.5 ^ A y = 1.0 
0.5 A y = 1 • o Sj 


One-dimensional Fields 



XX 


ZX 

xz. 
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h = altitude of bird 
S = skin depth 

a = first layer thickness for a two-layer model (not relevant 
in these calculations) 

T , t and T are generally referred to as the Sommerfeld 
0 ' 1 ' 2 

integrals and must be evaluated numerically. Thus far we have 
coded only T Q . 
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Note: The contribution of vertical surface elements to the 

models that are used in these calculations is small compared to 
the horizontal elements and so has not been included in the 
section presenting results. 

With the detector located at the origin of the 
coordinate system, the transmitter has Cartesian coordinates 
(Xg,0,0) and an element in the source surface has Cartesian 
coordinates (x,y,-h). R is the complex distance between a source 
element and the detector, and R' is the complex distance between 
a source element and the transmitter. 
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results 


in this section we discuss the results of the 
investigations in the context of the tacks set forth in the 

contract • 

Task l: The implementation of a MIM half space mver 

algorithm for AEM data. 

This task has been completed, and a printout of the 
code and a magnetic tape containing the code has been delivered 

to CRREL. 

Task 2 ! implementation of a MIM algorithm to calculate 
approximate AEM signatures of ice keels as functions of their 

width and depth. 

This task has been completed. MIM two-dimensional 
computer codes for rectangular and triangular sea ice keel models 
have also been delivered to CREEL in both printout and magnetic 
rape forms. We were initially tasked to produce algorithms and 
sample calculations of the MIM two-dimensional ZZ fields 
(horizontal coplanar coil configuration, . We have extended that 
task to include the XX (horizontal coaxial coils) and ZX fields. 
These calculations have been applied to two ice keel models. For 
the CRREL model based on ground truth data of an ice keel in 
Prudhoe Bay. the ZZ and XX fields have been inverted using the 
one-dimensional MIM inversion to produce values for the model 
parameters. The Berkeley model (Becker et al., 1987) is a 
triangular ice keel for which the XX field has been calculated 
using a finite element method. We compare results of these two 
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calculations. 


A detailed discussion follows. 


Task 3 : implementation of an analytic continuation algorithm of 

AEM data . 

The feasibility of analytically continuing AEM data up 
and down 10 m by a Taylor's series expansion of the MIM field has 
been demonstrated. A preliminary report of these results was 
made at the American Geophysical Union meeting in San Francisco 
in Dec 1988, and the abstract published (Bergeron et al., 1988). 
A more detailed discussion and illustration of this work follows 


below . 


Task 4: Deconvolution algorithm for ice keel signal signatures. 

The results for the preceding tasks show that the two- 
dimensional MIM ice keel fields (with the exception of the ZX and 


probably the XZ fields) track the one-dimensional fields in a 
"smoothed out" way. It is anticipated that an efficient 
deconvolution algorithm of the two-dimensional fields by a line 
impulse source signal will sharpen the two-dimensional fields and 
bring them into closer agreement with the one-dimensional fields, 
thereby bringing one-dimensional inversion results of the 
dimensional fields into closer agreement with the input model. 
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TASK TWO 


Approximate ZZ MIM Signal Signatures 
of Rectangular and Triangular Ice Keels 

The sea ice is assumed to be transparent to the AEM 
field, i.e., its conductivity is assumed negligible. All lengths 
in these models are scaled to the AEM skin depth in sea water and 
all fields calculated are for the horizontal coplanar ZZ coil 
configuration. 

Rectangular Keel 

bird 

t 

H 

sea ice I 

i S 

| DK 

1 1 

| seafwater 

< WK V 
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Definition of variables: 

H = altitude of bird above sea water/ ice interface 

WK = width of ice keel 
DK = depth of ice keel 

DX = numerical integration increment parallel to survey path 
D Y „ numerical integration increment perpendicular to survey path 
xo - coil spacing in AEM bird (for example, if the coil spacing 
is 6 m and XO is 0.6, then the shin depth is 10 m, the 
altitude of the bird above the seawater is 3.0 m if H = 
3.0, and the width and depth of the keel are 10 m if WK 

and DK are 1.0) 

in Figures 2-1 through 2-11, H Z /H P is plotted versus 
range (in skin depths) for rectangular and triangular keel 
models. The keels are centered at aero for various combinations 
Of WK, DK, and H. The solid and dashed curves are the rea 
quadrature signals, respectively. The larger pairs of signals 


1A 




are produced by the rectangular keels. 

in Figures 2-12 and 2-13, AMP is the maximum signal in 

ppm produced by the ice keel as the bird passes overhead. Figure 
2-12 shows AMP versus WK, and of course the signal saturates at 
about 600 ppm. This is the difference in the signals produced at 
altitudes of H = 3 and H+DK = 4. Figure 2-13 shows the variation 
in AMP with keel depth DK for a fixed keel width. The saturated 

value of AMP in this case is about 140 ppm. 

in Figures 2-14, 2-15, and 2-16 the 50% signal level 

signature width (WIDTH) is plotted versus keel width for constant 
keel depth, or versus keel depth for constant keel width. Note 
that the residual 50% signal width for narrow ice keels (WK < 4) 
is approximately 7 for H = 3. Thus narrow ice keels produce a 
50% width signal approximately equal to twice the bird altitude 

plus ice thickness. 
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Berkeley Triangular Keel Model 

The XX response of the triangular keel model shown in 
the bottom of Fig. B-l was calculated by Becker et ml. (1985) by 
means of a finite element algorithm. The results of the 
calculation at ten survey points are shown in the upper half of 
pig, B-l. This numerical procedure reportedly reguired 30 
minutes on a Cray supercomputer. The range of the calculation 
only approximately three times the keel width-not enough to 
reach the homogeneous halfspace values of the sommerfeld one- 
dimensional field of -312 ppm inphase, -62 quadrature. The same 
eodel is shown in Fig. B-2 but range, keel width, and keel depth 
are displayed in units of skin depth, which is 5.03 m given the 
assumed values for sea conductivity and bird transmitter 
frequency. The HIM one-dimensional and approximate two- 
dimensional fields for the XX coil configuration are shown in 
Fig. B-3 for the same model parameters. The first feature to 
note in these curves is the displacement of the minimum in the 
two-dimensional fields relative to the one-dimensional fields. 
This is caused by our computational scheme, and since this 
discrepancy continually reoccurs a brief explanation is required. 
The minimum values for the two-dimensional fields occur at about 
-6 m on the range scale. This is the approximate coordinate of 
the receiver coil when the transmitting coil is at the origin, 
which is the coordinate of the keel bottom relative to the model. 
At this location of the transmitter, the smallest currents are 
induced in the model and hence a minimum signal is detected at 
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the receiver coil. On the other hand, in the one-dimensional 
case, the altitude of the bird above the homogeneous half space is 
taken to be the vertical separation of the receiver coil above 
the model. This distance is a maximum when the receiver coil is 
directly above the bottom of the keel. Thus it is at this 
location where the one-dimensional field minimum occurs. 

The MIM XX two-dimensional field is in fair agreement 
with the Berkeley calculation. The signal widths of the two 
curves are approximately the same. The minimum field for the MIM 
two-dimensional calculation is intermediate between the Berkeley 
two-dimensional calculation and the MIM one-dimensional field 

minima. 

Figure B-4 shows the two-dimensional MIM fields for the 
ZZ coil configuration. The signal is broader and the maximum 
change in signal is smaller than for the XX configuration. Ho ZZ 

field was reported in the Berkeley report. 

Finally, the ZX MIM one- and two-dimensional fields are. 

displayed in Fig. B-5. This figure clearly shows that the MIM 
approximation predicts that the ZX field produces a Xeel signature 
that is very different from the one-dimensional signal. The 
leading edge downslope in the model produces a tracking downslope 
in the one-dimensional fields but a large positive hump two- 
dimensional signal and a nearly mirror negative hump indicative 
of the keel upslope. The peak to peak separation is 
approximately 25 m, which is comparable to the keel width of 18 

m. 
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CRREL Model 


Figure C-l through C-9 are the results of ground truth 
measurements of a Frudhoe Bay ice keel by CRREL personnel and 
contractors. The ordinant of the figure is the total ice/snow 
thickness, i.e., snow cover plus freeboard plus ice keel. A 
uniform thickness region has been added to each end of the keel 
in order to demonstrate that the two-dimensional MIM calculations 
of AEM fields are in agreement with the one-dimensional field in 
the uniform thickness regions of the model. Figure C-2 shows the 
altimeter reading of the bird-to-snow vertical distance measured 
during a helicopter traverse of the keel. The final model 
parameter used in these calculations is the sea conductivity, a 
= 3.1 S/m. 

Figures C-3, C-4, and C-5 show the results of one- and 
two-dimensional calculations for the ZZ coil configuration 
(horizontal coplanar) . Figure C-3 shows the inphase and 
quadrature one- and two-dimensional ZZ fields. Both fields track 
the model but the two-dimensional results are smoother and the 
variations are smaller than the one-dimensional results. Figures 
C-4 and C-5 show the results of the one-dimensional MIM inversion 
of the one- and two-dimensional fields. Figure C-4 gives the ice 
thickness (inverse distance from bird to sea surface minus the 
laser altimeter reading) and Figure C-5 the inversion results for 
the sea conductivity. The inversion results for the one- 
dimensional fields are in close agreement with the input model 
whereas the two-dimensional results indicate a smoother keel no 
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as dependent on the model. 

Figures C-6, C-7, and C-8 show results of the sane 

calculations for the XX configuration (horizontal coaxial) . A 

comparison of Figures C-4 and C-6 shows that the XX two- 

dimensional field to be closer to its one-dimensional 

• * f ha n is true for the two- and one-dimensional 

counterpoint than is true i 

_ „ This observation is confirmed 

fields for the ZZ configuration. This on 

results (ice thickness and conductivi y) ° 
bv the MIM inversion results n 

. o n and r-a Figure C-9 shows the 

the XX fields shown in Figures 

one - and two-dimensional fields for the mixed coil configurations 
Z X (vertical transmitter dipole/horizontal receiver coil axis) . 

Fne most striking feature of these curves is the large difference 
between the two- and one-dimensional fields at the leading e ge 
of the keel. The two-dimensional results show an enhanced 
secondary field whereas the one-dimensional field falls off with 
increasing distance between bird and sea surface. Thus a Z 
field would seem to provide a signature for a sudden increase in 
ice thickness. In this model the effect on the ZX field of the 
sudden decrease in ice thickness at about 100 m downrange (Figure 

C-!) tends to be canceled by the increasing altitude o 

, ... (Figure C-2) . Nevertheless 

over the same portion of the range (ng 

the two-dimensional fields do show a relatively sharp decrease 
from SO to 110 m: a mirror image trailing edge response to the 

initial leading edge response. This is more clearly seen in the 
ZX field of Figure C-10 in which the model has been modified to a 

constant bird altitude of 18 m. 
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in principle these departures of the two-dimensional 
calculation for the ZX configuration from the one-dimensional ZX 
fields can be compared to the measured ZX response of scaled edge 
models (Dallal, 1985). Dallal measured the ZX field in the tame 
domain over a sheet of brass as a function of distance from the 
edge of the sheet. It is possible for us to calculate the ZX 
field for such mode) s in the frequency domain over a sufficiently 
large range of frequencies and then Fourier transform these 
results for ZX(») into ZX(t), thereby allowing a comparison of 
the HIM approximation of ZX(t) with the scale model measurements. 
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TASK THREE 


Analytic Continuation of AEM Signal 
Figure 3-1 is a summary sheet of the MIM model along 
with the indicated Taylor's series expansion of the vertical 
component Z of the normalized secondary field (Bergeron et al., 
1988) . Figure 3-2 shows the smoothed laser altimeter versus 
range of a Prudhoe Bay survey line. It, along with an assumed 
conductivity of 2.7 S/m, constitutes a uniform half space model. 

The Sommerfeld integral expression, T 0 , for a secondary field is 
employed with this model to generate the inphase (real) and 
quadrature (imaginary) fields shown in Figs. 3-3 and 3-4. 

Figures 3-5 and 3-6 show the absolute percent difference between 
the Sommerfeld field calculated directly at a 40 m altitude and 
the fields shown in Fig. 3-3 and 3-4 and which are analytically 

continued to h = 40 m, i.e., 

% | AZ/Z | = | { [ Z s ( 40) - Z cont to 40 (b) ]/Zg(40) } X 100 

Figures 3-7 and 3-8 show the percent errors that result from an 
upward continuation of a signal from 30 m as a function of 
altitude for skin depths of 5 and 25 m, respectively. Figures 
3-9 and 3-10 show the corresponding percent errors that result 
from a downward continuation of a signal from 50 m. All of the 
figures indicate that a smaller error occurs for an upward 
continuation than a downward continuation. The utility of the 
technique depends on the subtlety of the anomaly that one is 
searching for in the data. Only if the error produced by the 
continuation of the data to a fixed reference altitude is less 
than the anticipated anomaly signal will this procedure be useful 
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TASK FOUR 


Deconvolution of a keel signal 


Figure 4—1 shows the one— and two-dimensional inphase 
and quadrature fields produced for a rectangular keel model of 
width 90 m and depth 5 m. Figure 4-2 shows the inphase one- and 
two-dimensional fields, again along with the result of a 
deconvolution of the two-dimensional field by a line impulse 
function. The signal width of the deconvolved fields is 
narrowed, which is the desired result. Gibbs oscillations are 
introduced by the deconvolution process. Similar results are 
shown in Figure 4-3 for the quadrature component of the fields. 
The step in the deconvolved two-dimensional quadrature field at 
about 20 m is probably an artifact of the line impulse function 
used in the calculation. 

Figures 4-4 through 4-7 show the results of a one- 
dimensional inversion of the fields displayed in Figures 4-1 
through 4—3. Figures 4—4 and 4—5 illustrate ice thickness and 
ice conductivity results for the one- and two-dimensional fields, 
and Figures 4—6 and 4—7 show the one— dimensional results, again 
along with those for the two-dimensional deconvolved fields. 

The deconvolved two-dimensional inversion results are 
not noticeably different from those of the original two- 
dimensional fields. The chief benefit of the deconvolution 
process that can be seen in this example is the narrowing of the 
keel signal. Further investigation with a variety of models 
should be pursued. 
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CONCLUSION 


The manner in which most of the tasks were completed 
exceeded initial specifications, e.g., the calculation of the XX 
and ZX ice keel fields. 

The validation of the MIM two- and three-dimensional 
fields by comparison with the results of more accurate (but more 
CPU time consuming) numerical calculations and model measurements 
should be the primary thrust in a continuing investigation. 
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Sea Ice Inversion 

The MIM inversion of sea ice AEM data taken at two 
frequencies, 1000 Hz (f lo ) and 250 kHz (f hi ) , proceeds as 
follows . 

The low frequency data are first inverted to give 
(h + dp) and a 2 , where h is the altitude of the bird above the sea 
ice, d is the ice thickness, and a 2 is the electrical 
conductivity of the sea. This inversion assumes the following: 

1. The skin depth of the sea ice Sp(lo) is much greater than the 
thickness of the sea ice and hence the sea ict; is effectively 
transparent to the low frequency primary signal. 

2. The sea bottom does not affect the secondary field. This 
assumption is valid provided that the sea depth d 2 is greater 
than twice the low frequency skin depth of the sea, i.e. , d 2 > 2 
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52(1°) • These two assumptions allow a halfspace inversion. The 
algebra and computer algorithms for the halfspace inversion are 
given in the Appendix. 

It is assumed that the altitude h is independently 
determined by a radar or laser altimeter. Thus the inversion 
results in a local value for the sea ice thickness d^ and the 
conductivity of the sea water, <72- These results are employed in 
the inversion of the high frequency data to determine the sea ice 
conductivity . 


Outline of high frequency inversion 

First a half space inversion of the high frequency data 
is performed. This produces an effective skin depth S e ff which 
lies in the range 5 2 < S e ff < $]_, and is a function of ice 
thickness d^. The effective high frequency skin depth is 
combined with the altimeter reading h to form the ratio A eff - 
2h/6 e ff . The ad hoc normalization function employed in MIM 
inversion is a function of A e ff, i.e., 

Zj^im “ Z s (normalized) — F(A e ££) 

For d-L « 6^, $ e ff ~ & 2 an< * ^ or **1 > ^1’ t * ien ^eff - 
6-^. Since d^ is known from the low frequency Inversion, this 
latter case may be recognized and hence the first layer 
conductivity o-y is determined from 6 e ff by ^ 

o 1 - 2/[hq f (hi) S eff 2 ] , V 

where /iq is the vacuum magnetic permittivity. The condition 
di > 2 S-l only occurs 'or thick (dj^ > 10m) , highly conducting (c^ 
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> 0.027 S/m) sea ice. For the intermediate, more general 

situation where 5 e ff/ d l is of the ° rder ° f Unity ’ inversi ° n 

procedure to be used is that described below. 

The MIM relationship between the complex two-layer 

correction factor Q and the high frequency AEM field is 
algebraically transformed into two simultaneous transcendental 
real equations with argument d^, where 6 1 is the unknown 
quantity. All other quantities in these equations are known. 

Each of these equations has in general several roots, BUT only one 
common root. The explicit functions that occur respectively in 
these equations are tan(d l/5l ) and tanh^/S,) . A root-finding 
algorithm is first applied to the tan^/^) equation. When a 
root is determined, that root is inserted in the tanh(d 1 /« 1 ) 
equation to test if it is also a root of the tanh(d 1 /5 1 ) 
equation. If not, the algorithm continues in its determination 
of the real roots of the tan (d^) equation until the root is 
found that simultaneously satisfies both equations. The first 
layer skin depth 6 ± and conductivity a x are given by that 
simultaneous root. 

The range of applicability of the root finding 
algorithm is given by 0 02 < dp/Jp <2.5. These limits can be 
understood in physical terms. For dp/6p > 2.5 the sea ice is 
effectively a halfspact as has been already noted, and a two 
layer model is inappropriate. For dp/Sp < 0.02 the perturbation 
produced on the secondary AEM field by the sea ice cover is lost 
in the computer -noise- caused by roundoff, etc., and will 
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certainly be undetectable in the noise and drift present in even 
ideal real data, where noise and drift are greater than about 1 

ppm. 

The lower ice. thickness limit on the detectability of 
sea ice conductivity is illustrated in the following table which 
assumes a value of sea water conductivity of o 2 " 27 S/m and an 
operating frequency of 250 Khz. 

5^ minimum d-^ 

-6m - 0.1 m 

- 20 m - 0.3 m 

The algebraic details of this procedure and the root finding 
algorithm are given in the Appendix. 

Results 

The MIM inversion procedure that has been described is 
applied to several sea ice models. In all of the models used the 
low and high frequencies assumed for the AEM system are 1 kHz and 
250 kHz, respectively; the altitude of the AEM bird is 25 m; the 
conductivity of the sea water a 2 is 2.7 S/m; and the conductivity 
of the sea ice a x for each model has input values of 0.027 S/m, 
0.0054 S/m, and 0.0027 S/m. Thus the ratio K of the 
conductivities of sea water to sea ice has the values 100, 500, 

and 1000, respectively. 

With these general conditions, the first model of ice 


a 2 /^i 

100 

1000 


4 



thickness versus range (fiducial number) is given in Figure 1. 

The ice thickness increases linearly with increasing range. The 
results of the inversion for o ^ are shown in Figure 2. The 
inversion values for a ] are in fair agreement with the input 
values except for the case with a ^ — 0.027 S/m. The problem 
occurs at an ice thickness of approximately 9.5m. For a ^ - 
0.027 S/m the skin depth of the sea ice is about 6 m, thus the 
ratio of ice thickness to skin depth (which is the argument of 
both the tan and tanh junctions) is about pi/2, where the tangent 
becomes singular and double valued. More importantly, in the 
immediate vicinity of pi/2, tan(d \/S^) varies rapidly. In spite 
of this, the root finding inversion algorithm still works when 
the exact forward MIM field ZZ(MIM) is used as the input field 
(see Table 1). When a simultaneous root cannot be found for the 
normalized Sommerfeld field in the vicinity of n/2 , a value of 
1.55 is assumed for x. See Table 2. It is the residual 
differences between the normalized Sommerfeld field (or real 
field data) and the exact MIM field that causes the root finding 
algorithm that we are currently using to fail for x - d^/S^ “ 
ar/2. 

It should be noted that this value* of n/2 will most 
likely not be encountered in field surveys where ice 
conductivities will generally be less than 0.0054 S/m (K - 500). 
Table 3 shows that for K - 500, x is less than n/2 for sea ice 
thicknesses up to 20 m. 

A shallow ice keel model is shown in Figure 3. Figure 4 
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shows the values of o x for this model produced by the inversion 
algorithm for K - 100, 500, and 1000. These results are also 

listed in Tables 5, 6, and 7. 

In all of the tables we have included the results of 

the halfspace inversion of the high frequency data which gives 

a eff* Xt can be S6en fcr the CaSe K “ 100 When X > 2A ’ ° eff ~ 

This demonstrates that when the ice thickness is greater 

° input • 

than 2.4 skin depths, a halfspace inversion yields good results 
for the ice conductivity. Although values of x > 2.4 will 
probably not be found in survey data taken at a high frequency of 
250 kHz, still higher frequencies of about 1 MHz will bring x 
into this range. 

Finally, the results of the inversion are shown in 
Figure 3 for a shallow ice keel model. The tabulated results are 

shown in Tables 5, 6, and 7. 

In summary, the present inversion algorithm for a 1 

works well except in the vicinity of dj/^ - */2. We are 
continuing efforts to modify and improve the existing algorithm. 
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SIGMA1 

(ACTUAL) =0.027 S/M 


WORKING WITH 



ZZMIM 





XSTART= 

= 0.5000000 



KNOW D1 AND 

SIG2, SOLVE 

FOR X=D1/DELT1 AND HENCE SIG1. 

SIGeff 

IS THE HALF -SPACE 

: EFFECTIVE 

CONDUCTIVITY, 

AND IT 

IS COMPUTED ONLY 

WHEN INVERTING ZZSOM 

FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0156 

0.0247 

2.6960 

1.0 

0.5 

0.0781 

0.0247 

2.5941 

2.0 

1.0 

0.1563 

0.0247 

2.2813 

3.0 

1.5 

0.2344 

0.0247 

1.7051 

4.0 

2.0 

0.3250 

0.0268 

1.1303 

5.0 

2.5 

0.4063 

0.0268 

0.7022 

6.0 

3.0 

0.4875 

0.0268 

0.4009 

7.0 

3.5 

0.5687 

0.0268 

0.2198 

8.0 

4.0 

0.6500 

0.0268 

0.1332 

9.0 

4.5 

0.7312 

0.0268 

0.0884 

10.0 

5.0 

0.8125 

0.0268 

0.0613 

11.0 

5.5 

0.9000 

0.0271 

0.0461 

12.0 

6.0 

0.9818 

0.0271 

0.0373 

13.0 

6.5 

1.0636 

0.0271 

0.0315 

14.0 

7.0 

1.1455 

0.0271 

0.0277 

15.0 

7.5 

1.2273 

0.0271 

0.0253 

16.0 

8.0 

1.3091 

0.0271 

0.0238 

17.0 

8.5 

1.3909 

0.0271 

0.0228 

18.0 

9.0 

1.4727 

0.0271 

0.0224 

19.0 

9.5 

1.5500 

0.0270 

0.0221 

20.0 

10.0 

1.6355 

0.0271 

0.0222 

21.0 

10.5 

1.7172 

0.0271 

0.0223 

22.0 

11.0 

1.7990 

0.0271 

0.0226 

23.0 

11.5 

1.8808 

0.0271 

0.0230 

24.0 

12.0 

1.9626 

0.0271 

0.0234 

25.0 

12.5 

2.0443 

0.0271 

0.0237 

26.0 

13.0 

2.1218 

0.0270 

0.0241 

27.0 

13.5 

2.2035 

0.0270 

0.0244 

28.0 

14.0 

2.2851 

0.0270 

0.0247 

29.0 

14.5 

2.3667 

0.0270 

0.0250 

30.0 

15.0 

2.4483 

0.0270 

0.0252 

31.0 

15.5 

2.5299 

0.0270 

0.0254 

32.0 

16.0 

2.6115 

0.0270 

0.0255 

33.0 

16.5 

2.6931 

0.0270 

0.0257 

34.0 

17.0 

2.7747 

0.0270 

0.0258 

35.0 

17.5 

2.8563 

0.0270 

0.0258 

36.0 

18.0 

2.9379 

0.0270 

0.0259 

37.0 

18.5 

3.0195 

0.0270 

0.0259 

38.0 

19.0 

3.1012 

0.0270 

0.0259 

39.0 

19.5 

3.1828 

0.0270 

0.0259 

40.0 

20.0 

3.2644 

0.0270 

0.0259 




TABLE 2 


SIGMA1 (ACTUAL) =0.027 S/M 
WORKING WITH 


ZZSOM 

ASi 5 s?G°!°S0LVE FOR X-D1/DELT1 AND HENCE 8101. 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVITY, 
and it is computed only when INVERTING ZZSOM 


FID 

D1 X 

0.0 

0.1 0.0156 

1.0 

0.5 0.0781 

2.0 

1.0 0.1563 

3.0 

1.5 0.2344 

4.0 

2.0 0.3375 

5.0 

2.5 0.4087 

6.0 

3.0 0.4904 

7.0 

3.5 0.5722 

8.0 

4.0 0.6664 

9.0 

4.5 0.7497 

10.0 

5.0 0.8330 

11.0 

5.5 0.9163 

12.0 

6.0 0.9996 

13.0 

6.5 1.0829 

14.0 

7.0 1.1662 

15.0 

7.5 1.2495 

16.0 

8.0 1.5500 

17.0 

8.5 1.5500 

18.0 

9.0 1.5500 

19.0 

9.5 1.4875 

20.0 

10.0 1.5500 

21.0 

10.5 1.7131 

22.0 

11.0 1.7947 

23.0 

11.5 1.8763 

24.0 

12.0 1.9579 

25.0 

12.5 2.0394 

26.0 

13.0 2.1210 

27.0 

13.5 2.2026 

28.0 

14.0 2.2842 

29.0 

14.5 2.3657 

30.0 

15.0 2.4473 

31.0 

15.5 2.5289 

32.0 

16.0 2.6105 

33.0 

16.5 2.6921 

34.0 

17.0 2.7736 

35.0 

17.5 2.8552 

36.0 

18.0 2.9383 

37.0 

18.5 3.0199 

38.0 

19.0 3.1024 

39.0 

19.5 1.5500 

40.0 

20.0 1.5500 


SIG1 

SIGeff 

0.0247 

2.6973 

0.0247 

2.5953 

0.0247 

2.2816 

0.0247 

1.7033 

0.0289 

1.1265 

0.0271 

0.6976 

0.0271 

0.3973 

0.0271 

0.2178 

0.0281 

0.1324 

0.0281 

0.0884 

0.0281 

0.0619 

0.0281 

0.0471 

0.0281 

0.0386 

0.0281 

0.0330 

0.0281 

0.0294 

0.0281 

0.0270 

0.0380 

0.0256 

0.0337 

0.0247 

0.0301 

0.0242 

0.0248 

0.0240 

0.0243 

0.0240 

0.0270 

0.0241 

0.0270 

0.0244 

0.0270 

0.0247 

0.0270 

0.0250 

0.0270 

0.0253 

0.0270 

0.0256 

0.0270 

0.0259 

0.0270 

0.0261 

0.0270 

0.0264 

0.0270 

0.0265 

0.0270 

0.0267 

0.0270 

0.0268 

0.0270 

0.0269 

0.0270 

0.0270 

0.0270 

0.0271 

0.0270 

0.0271 

0.0270 

0.0271 

0.0270 

0.0271 

0.0064 

0.0271 

0.0061 

0.0271 




TABLE 3 


SIGMA1 (ACTUAL) =0.0 054 S/M 

WORKING WITH 

ZZSOM 

XSTART= 0.5000000 

KNOW D1 AND SIG2 , SOLVE FOR X=D1/DELT1 AND HENCE SIG1. 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVITY, 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0078 

0.0062 

2.6998 

1.0 

0.5 

0.0391 

0.0062 

2.6784 

2.0 

1.0 

0.0781 

0.0062 

2.6794 

3.0 

1.5 

0.1172 

0.0062 

2.5616 

4.0 

2.0 

0.1563 

0.0062 

2.3817 

5.0 

2.5 

0.1953 

0.0062 

2.1671 

6.0 

3.0 

0.2344 

0.0062 

1.7031 

7.0 

3.5 

0.2734 

0.0062 

1.1052 

8.0 

4.0 

0.2930 

0.0054 

0.7481 

9.0 

4.5 

0.3296 

0.0054 

0.5329 

10.0 

5.0 

0.3662 

0.0054 

0.3495 

11.0 

5.5 

0.4028 

0.0054 

0.2399 

12.0 

6.0 

0.4395 

0.0054 

0.1807 

13.0 

6.5 

0.4886 

0.0057 

0.1303 

14.0 

7.0 

0.5262 

0.0057 

0.0979 

15.0 

7.5 

0.5637 

0.0057 

0.0723 

16.0 

8.0 

0.6076 

0.0058 

0.0536 

17.0 

8.5 

0.6455 

0.0058 

0.0414 

18.0 

9.0 

0.6835 

0.0058 

0.0325 

19.0 

9.5 

0.7215 

0.0058 

0.0264 

20.0 

10.0 

0.7595 

0.0058 

0.0213 

21.0 

10.5 

0.7974 

0.0058 

0.0178 

22.0 

11.0 

0.8354 

0.0058 

0.0148 

23.0 

11.5 

0.8734 

0.0058 

0.0126 

24.0 

12.0 

0.9114 

0.0058 

0.0109 

25.0 

12.5 

0.9493 

0.0058 

0.0097 

26.0 

13.0 

0.9873 

0.0058 

0.0087 

27.0 

13.5 

1.1253 

0.0070 

0.0082 

28.0 

14.0 

1.1670 

0.0070 

0.0075 

29.0 

14.5 

1.2086 

0.0070 

0.0070 

30.0 

15.0 

1.2503 

0.0070 

0.0066 

31.0 

15.5 

1.2112 

0.0062 

0.0063 

32.0 

16.0 

1.2503 

0.0062 

0.0060 

33.0 

16.5 

1.2894 

0.0062 

0.0058 

34.0 

17.0 

1.3285 

0.0062 

0.0056 

35.0 

17.5 

1.3675 

0.0062 

0.0055 

36.0 

18.0 

1.4066 

0.0062 

0.0054 

37.0 

18.5 

1.5500 

0.0071 

0.0053 

38.0 

19.0 

1.5500 

0.0067 

0.0052 

39.0 

19.5 

1.5500 

0.0064 

0.0052 

40.0 

20.0 

1.5500 

0.0061 

0.0051 




TABLE 4 

f 




SIGMA1 

(ACTUAL) =0.0027 S/M 


WORKING WITH 



ZZSOM 





XSTART= 

= 0.5000000 



KNOW D1 AND 

SIG2. SOLVE 

FOR X=D1/DELT1 AND HENCE SIG1. 

SIGeff 

IS THE HALF -SPACE 

: EFFECTIVE 

CONDUCTIVITY , 

AND IT 

IS COMPUTED ONLY 

WHEN INVERTING ZZSOM 

FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0078 

0.0062 

2.7001 

1.0 

0.5 

0.0391 

0.0062 

2.6894 

2.0 

1.0 

0.0781 

0.0062 

2.7387 

3.0 

1.5 

0.0586 

0.0015 

2.7221 

4.0 

2.0 

0.0781 

0.0015 

2.7030 

5.0 

2.5 

0.0977 

0.0015 

2.7093 

6.0 

3.0 

0.1672 

0.0031 

2.3818 

7.0 

3.5 

0.1951 

0.0031 

1.6995 

8.0 

4.0 

0.2229 

0.0031 

1.2689 

9.0 

4.5 

0.2508 

0.0031 

1.0030 

10.0 

5.0 

0.2786 

0.0031 

0.7038 

11.0 

5.5 

0.2874 

0.0028 

0.5150 

12.0 

6.0 

0.3135 

0.0028 

0.4217 

13.0 

6.5 

0.3396 

0.0028 

0.3189 

14.0 

7.0 

0.3657 

0.0028 

0.2518 

15.0 

7.5 

0.3918 

0.0028 

0.1897 

16.0 

8.0 

0.4180 

0.0028 

0.1405 

17.0 

8.5 

0.4566 

0.0029 

0.1086 

18.0 

9.0 

0.4835 

0.0029 

0.0843 

19.0 

9.5 

0.5103 

0.0029 

0.0680 

20.0 

10.0 

0.5372 

0.0029 

0.0528 

21.0 

10.5 

0.5640 

0.0029 

0.0426 

22.0 

11.0 

0.5909 

0.0029 

0.0332 

23.0 

11.5 

0.6177 

0.0029 

0.0265 

24.0 

12.0 

0.6446 

0.0029 

0.0217 

25.0 

12.5 

0.6715 

0.0029 

0.0180 

26.0 

13.0 

0.6983 

0.0029 

0.0151 

27.0 

13.5 

0.7252 

0.0029 

0.0128 

28.0 

14.0 

0.7520 

0.0029 

0.0109 

29.0 

14.5 

0.7789 

0.0029 

0.0095 

30.0 

15.0 

0.8058 

0.0029 

0.0083 

31.0 

15.5 

0.9326 

0.0037 

0.0073 

32.0 

16.0 

0.9627 

0.0037 

0.0067 

33.0 

16.5 

0.9928 

0.0037 

0.0061 

34.0 

17.0 

1.0229 

0.0037 

0.0055 

35.0 

17.5 

1.0529 

0.0037 

0.0051 

36.0 

18.0 

1.0830 

0.0037 

0.0047 

37.0 

18.5 

1.1131 

0.0037 

0.0044 

38.0 

19.0 

1.1432 

0.0037 

0.0041 

39.0 

19.5 

1.1733 

0.0037 

0.0038 

40.0 

20.0 

1.2034 

0.0037 

0.0036 




TABLE 5 


=SSS — = — — 

— 


SIGMA1 

(ACTUAL) =0.027 S/M 


WORKING 
ZZSOM 
XSTART= 
KNOW D1 

WITH 

0.5000000 
AND SIG2 , SOLVE 

FOR X=D1/DELT1 AND HENCE SIG1. 

SIGeff 

IS THE HALF -SPACE 

i EFFECTIVE 

CONDUCTIVITY, 

AND IT 

IS COMPUTED ONLY 

WHEN INVERTING ZZSOM 

FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0156 

0.0247 

2.6973 

1.0 

0.5 

0.0781 

0.0247 

2.5953 

2.0 

1.0 

0.1563 

0.0247 

2.2046 

3.0 

1.5 

0.2469 

0.0274 

1.6041 

4.0 

2.0 

0.3292 

0.0274 

1.0232 

5.0 

2.5 

0.4115 

0.0274 

0.6057 

6.0 

3.0 

0.4938 

0.0274 

0.3527 

7.0 

3.5 

0.5760 

0.0274 

0.2107 

8.0 

4.0 

0.6583 

0.0274 

0.1324 

9.0 

4.5 

0.7469 

0.0279 

0.0884 

10.0 

5.0 

0.8299 

0.0279 

0.0631 

11.0 

5.5 

0.9191 

0.0283 

0.0481 

12.0 

6.0 

1.0027 

0.0283 

0.0388 

13.0 

6.5 

1.0862 

0.0283 

0.0330 

14.0 

7.0 

1.1652 

0.0281 

0.0294 

15.0 

7.5 

1.2484 

0.0281 

0.0270 

16.0 

8.0 

1.5500 

0.0380 

0.0256 

17.0 

8.5 

1.5500 

0.0337 

0.0247 

18.0 

9.0 

1.5500 

0.0301 

0.0242 

19.0 

9.5 

1.5000 

0.0253 

0.0240 

20.0 

LO.O 

1.6274 

0.0268 

0.0241 

21.0 

9.5 

1.5500 

0.0270 

0.0240 

22.0 

9.0 

1.5500 

0.0301 

0.0242 

23.0 

8.5 

1.5500 

0.0337 

0.0247 

24.0 

8.0 

1.5500 

0.0380 

0.0256 

25.0 

7.5 

1.2500 

0.0281 

0.0270 

26.0 

7.0 

1.1667 

0.0281 

0.0294 

27.0 

6.5 

1.0833 

0.0281 

0.0330 

28.0 

6.0 

1.0000 

0.0281 

0.0388 

29.0 

5.5 

0.9167 

0.0281 

0.0481 

30.0 

5.0 

0.8333 

0.0281 

0.0631 

31.0 

4.5 

0.7500 

0.0281 

0.0884 

32.0 

4.0 

0.6667 

0.0281 

0.1324 

33.0 

3.5 

0.5833 

0.0281 

0.2107 

34.0 

3.0 

0.5000 

0.0281 

0.3527 

35.0 

2.5 

0.4167 

0.0281 

0.6057 

36.0 

2.0 

0.3333 

0.0281 

1.0232 

37.0 

1.5 

0.2500 

0.0281 

1.6041 

38.0 

1.0 

0.1667 

0.0281 

2.2046 

39.0 

0.5 

0.0833 

0.0281 

2.5953 

40.0 

0.1 

0.016V 

0.0281 

2.6973 




TABLE 6 

===== 

======='== = = = 

== 


SIGMA1 

(ACTUAL) =0.0 054 S/M 


WORKING 
ZZSOM 
XSTART= 
KNOW D1 

WITH 

0.5000000 
AND SIG2, SOLVE 

FOR X=D1/DELT1 AND HENCE SIG1. 

SIGef f 

IS THE HALF-SPACE EFFECTIVE 

CONDUCTIVITY, 

AND IT 

IS COMPUTED ONLY 

WHEN INVERTING ZZSOM 

FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0078 

0.0062 

2.6998 

1.0 

0.5 

0.0391 

0.0062 

2.6784 

2.0 

1.0 

0.0781 

0.0062 

2.5848 

3.0 

1.5 

0.1172 

0.0062 

2.3914 

4.0 

2.0 

0.1563 

0.0062 

2.0990 

5.0 

2.5 

0.1953 

0.0062 

1.7409 

6.0 

3.0 

0.2344 

0.0062 

1.3686 

7.0 

3.5 

0.2734 

0.0062 

1.0287 

8.0 

4.0 

0.2930 

0.0054 

0.7481 

9.0 

4.5 

0.3296 

0.0054 

0.5329 

10.0 

5.0 

0.3662 

0.0054 

0.3759 

11.0 

5.5 

0.4028 

0.0054 

0.2650 

12.0 

6.0 

0.4395 

0.0054 

0.1879 

13.0 

6.5 

0.4886 

0.0057 

0.1348 

14.0 

7.0 

0.5262 

0.0057 

0.0979 

15.0 

7.5 

0.5637 

0.0057 

0.0723 

16.0 

8.0 

0.6076 

0.0058 

0.0544 

17.0 

8.5 

0.6455 

0.0058 

0.0417 

18.0 

9.0 

0.6835 

0.0058 

0.0325 

19.0 

9.5 

0.7215 

0.0058 

0.0259 

20.0 

10.0 

0.7595 

0.0058 

0.0210 

21.0 

9.5 

0.7215 

0.0058 

0.0259 

22.0 

9.0 

0.6835 

0.0058 

0.0325 

23.0 

8.5 

0.6500 

0.0059 

0.0417 

24.0 

8.0 

0.6118 

0.0059 

0.0544 

25.0 

7.5 

0.5646 

0.0057 

0.0723 

26.0 

7.0 

0.5269 

0.0057 

0.0979 

27.0 

6.5 

0.4893 

0.0057 

0.1348 

28.0 

6.0 

0.4517 

0.0057 

0.1879 

29.0 

5.5 

0.4140 

0.0057 

0.2650 

30.0 

5.0 

0.3764 

0.0057 

0.3759 

31.0 

4.5 

0.3387 

0.0057 

0.5329 

32 . 0 

4.0 

0.3011 

0.0057 

0.7481 

33.0 

3.5 

0.2635 

0.0057 

1.0287 

34.0 

3.0 

0.2253 

0.0057 

1.3686 

35.0 

2.5 

0.1882 

0.0057 

1.7409 

36.0 

2.0 

0.1506 

0.0057 

2.0990 

37.0 

1.5 

0.1129 

0.0057 

2.3914 

38.0 

1.0 

0.0753 

0.0057 

2.5848 

39.0 

0.5 

0.0376 

0.0057 

2.6784 

40.0 

0.1 

0.0075 

0.0057 

2.6998 




TABLE 7 


::=r===;== 

=== 


SIGMA1 

(ACTUAL) =0.0027 S/M 


WORKING 
ZZSOM 
XSTART= 
KNOW D1 

WITH 

0.5000000 
AND SIG2. SOLVE 

FOR X=D1/DELT1 AND HENCE SIG1. 

SIGeff 

IS THE HALF-SPACE EFFECTIVE 

CONDUCTIVITY , 

AND IT 

IS COMPUTED ONLY 

WHEN INVERTING ZZSOM 

FID 

D1 

X 

SIG1 

SIGeff 

0.0 

0.1 

0.0078 

0.0062 

2.7001 

1.0 

0.5 

0.0391 

0.0062 

2.6894 

2.0 

1.0 

0.0781 

0.0062 

2.6414 

3.0 

1.5 

0.0586 

0.0015 

2.5379 

4.0 

2.0 

0.0781 

0.0015 

2.3690 

5.0 

2.5 

0.0977 

0.0015 

2.1376 

6.0 

3.0 

0.1672 

0.0031 

1.8600 

7.0 

3.5 

0.1951 

0.0031 

1.5616 

8.0 

4.0 

0.2229 

0.0031 

1.2689 

9.0 

4.5 

0.2508 

0.0031 

1.0030 

10.0 

5.0 

0.2786 

0.0031 

0.7760 

11.0 

5.5 

0.2874 

0.0028 

0.5912 

12.0 

6.0 

0.3135 

0.0028 

0.4462 

13.0 

6.5 

0.3396 

0.0028 

0.3353 

14.0 

7.0 

0.3657 

0.0028 

0.2518 

15.0 

7.5 

0.3918 

0.0028 

0.1897 

16.0 

8.0 

0.4180 

0.0028 

0.1437 

17.0 

8.5 

0.4566 

0.0029 

0.1096 

18.0 

9.0 

0.4835 

0.0029 

0.0843 

19.0 

9.5 

0.5103 

0.0029 

0.0654 

20.0 

LO.O 

0.5372 

0.0029 

0.0513 

21.0 

9.5 

0.5103 

0.0029 

0.0654 

22.0 

9.0 

0.4835 

0.0029 

0.0843 

23.0 

8.5 

0.4566 

0.0029 

0.1096 

24.0 

8.0 

0.4297 

0.0029 

0.1437 

25.0 

7.5 

0.4029 

0.0029 

0.1897 

26.0 

7.0 

0.3760 

0.0029 

0.2518 

27.0 

6.5 

0.3492 

0.0029 

0.3353 

28.0 

6.0 

0.3223 

0.0029 

0.4462 

29.0 

5.5 

0.2954 

0.0029 

0.5912 

30.0 

5.0 

0.2686 

0.0029 

0.7760 

31.0 

4.5 

0.2417 

0.0029 

1.0030 

32.0 

4.0 

0.2149 

0.0029 

1.2689 

33.0 

3.5 

0.1880 

0.0029 

1.5616 

34.0 

3.0 

0.1612 

0.0029 

1.8600 

35.0 

2.5 

0.1343 

0.0029 

2.1376 

36.0 

2.0 

0.1074 

0.0029 

2.3690 

37.0 

1.5 

0.0806 

0.0029 

2.5379 

38.0 

1.0 

0.0537 

0.0029 

2.6414 

39.0 

0.5 

0.0269 

0.0029 

2.6894 

40.0 

0.1 

0.0054 

0.0029 

2.7001 




Appendix 


The MIM representation of the normalized secondary 
field produced by induced ohmic currents in a two -layered 
conducting model (see Figure 1) for a horizontal coplanar coil 
pair is given by 

(H s /Hp) = ZZ - [2R 2 - 1]/[R 2 + 1] 5/2 (A-l) 

and 

R - [2h + (1-i) Q 8 x ]/p , (A-2) 

where [2h + (1-i) Q 6^] is the complex vertical distance 
separating the primary dipole source from the image source, h is 
the real altitude of the bird above the first layer surface, 

[(1-i) Q $]J/2 is the complex distance below the first layer 
surface of the image plane, p is the coil spacing, 5-^ the 
first layer (sea ice) skin depth, and finally Q is the two-layer 
correction factor given by 

q _ [{8^/62) + tanh{ ( l+i)d 1 /5 1 ) ]/[l + ( Sj /^) tanh{ (l+ijdj/^) ] . (A-3) 

d^ is the first layer thickness and 82 is the skin depth of the 
second layer (sea water). The thickness or depth of the second 
layer, & 2 » assumed to be greater than 2 6 2 *- n this analysis. 

(In order to determine the sea depth, a third lower frequency 
signal must be employed.) The R function can be interpreted 
geometrically as cotan <f>, where <f> is the complex angle indicated 
in Figure 1. 

In all MIM inversion schemes, Equation (A-2) is 
inverted by means of a polynomial expression which gives R as a 
function of ZZ: 
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R - 1/s - as - bs 1 - cs 5 - ds 7 , (A-4) 

where s - (ZZ/2) 1 / 3 , and for this coil configuration, a - 1, 
b - 9/8, c - 31/12, and d - 267/384. This inverted relationship 
is of paramount importance in all MIM inversion routines, i.e., 
from the value of R calculated from the AEM fields the values of 
the model parameters are determined. 

Halfspace inversion 

A halfspace inversion is defined by the condition Q - 
1. Thus Equation (A- 2) can be inverted to give values for the 
bird altitude h and the skin depth 

2h /p - R x + R 2 
S\/p " - ^2 

where R^ and R 2 are the real and imaginary components of R. In 
the event that the conditions which make Q — 1 are not satisfied 
(Q - 1 if ff 2 - or d l > 2 5^ , then the halfspace inversion 
results in an effective ;kin depth and bird altitude. 

Two layer inversion for sea ice conductivity 

In this inversion it is assumed that the only unknown 
is the first layer skin depth, i.e., the altitude, the second 
layer skin depth, and the first layer thickness are known from 
the altimeter reading and the low frequency inversion results. 

The algebraic manipulations of Equations (A- 2) and (A- 3) that 
result in the two simult aneous equations mentioned in the 
narrative proceed as follows: 
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First Q5^ is solved for explicitly from Equation (A- 2) 

to produce 

Q6 l - [p R(hi) - 2h]/(l - i) . (A-5) 

Next the expression for Q from Equation (A- 3) is substituted into 
Equation (A-5) and the resultant equation is rearranged to solve 
explicitly for tanh( (l+i)d^/6^) : 
tanh{ (l+i)d^/5^) - [ (C - 1) ( 5 2 / ^ l ) ] / [ 1 * C^/^l-l 

- D(6 1 ,6 2 c ) • (A ' 6) 

C is a knovm complex number given by 

C - [p R(hi) - 2 h]/[ (1 - i) 6 ? ] 
and D is a complex function of the unknown 6 ^ and the known 
quantities C and S 2 ■ If we set D - D x + i D 2 and expand 
tanhf (l+i)d^/5^) by 

tanh{(l+i)d 1 /5 1 ) - [tanh(di/«i) + i 6 ] / [\ + 

i tanh(d^/5^) tan(d^/6^)] 

we get 

[tanhCdj/S^ + i tan^/S^ ]/[l + i tarihCdj/fi].) tanCdj/fij.)] 

- + i D 2 ( a * 7 > 

Finally, if we equate the real and imaginary parts of Equation 

(A- 7) , we find 

tanhCdj/^) - D x - D 2 tanCd]/^) tariMd]/^) 
tanCd-L/S^ - D 2 + D x tanCd]/^) tanhCdj/S].) 

These can be combined to give 

D 2 tan^dj/fij.) - (D t 2 + D 2 2 - 1) tan(d 1 /5 1 ) - I) 2 - 0 (A-8a) 

and 

D x tanh 2 (d 1 /6 1 ) - (D x 2 + D 2 2 + 1) tanh^/S^ + D x - 0 . (A-8b) 
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2 

Both of these equations are of the quadratic form ax + b x + c, 
and hence explicit expressions for tanh(d 1 /6 1 ) and tan(d can 


be written. But first it should be noted that for the tan 
equation, a - - c - D 2 , and thus 

tanCdj/fi!) - - b/2a ± [(b/2a) 2 + 1] 1/2 , 

and similarly for the tanh equation, since a - c - D 1( 
tanh(d 1 /5 1 ) •» - b/2a + [(b/2a) 2 - 1] 1/2 

The physical constraint that the roots of these 
equations be real and positive results in the following 


solutions : 

tan(d 1 /6 1 ) =■ F 1 (d 1 /6 1 ) + [F 1 2 (d 1 /fi 1 ) + 1] 1/2 
for 0 < (d-L/5 < w/2, and 
tan(d !/«]_) - F 1 - [F x 2 + 1] 1/2 
for tt/ 2 < (d^/5]_) < where 

F x - [Dj 2 + D 2 2 - 1 ] / ( 2 D 2 ) 

tanh(d^/5^) - F 2 + [F 2 2 - 1]^^ 2 . 


(A-9a) 


(A- 9b) 


(A-10) 


where 

F 2 - [Dj^ +• D 2 2 + l]/(2 D x ) 

A root finding algorithm given below is used to find 
the real positive roots of Equation (A-9a) or (A-9b). These 
roots are substituted into Equation (A-10) to find the one root 
of Equation (A- 9a) or (A- 9b) that is simultaneously a root of 
Equation (A-10). This value of (d 1 /fi 1 ) is used to determine 
and in turn a 


10 




FIGURE 1 - Appendix 


RECEIVER TRANSMITTER 



* 


IMAGE 




Updated Status Report 


15 Mar 1990 
NASA Grant NAG- 1-804 

Determination of Design and Operation Parameters for Upper 
Atmospheric Research Instrumentation to Yield Optimum Resolution 

with Deconvolution 

Analysis of Response of the Arctic Sea Ice Environment 

to AEM Fields 

Clyde J. Bergeron and Juliette W. Ioup 
Department of Physics, University of New Orleans 

Sea Ice Inversion 

The MIM inversion of sea ice AEM data taken at two 
frequencies, 1000 Hz (f-^ 0 ) an< ^ 250 kHz (fj-ii^ > and for two AEM 
coil configurations, horizontal coplanar and vertical coaxial, 
is outlined below. 

The low frequency data are first inverted to give 
(h + d-j^ and o 2 , where h Ls the altitude of the bird above the sea 
ice, d is the ice thickness, and <? 2 is the electrical 
conductivity of the sea. This inversion assumes the following: 

1. The skin depth of the sea ice fi-^(lo) is much greater than the 
thickness of the sea ice and hence the sea ice is effectively 
transparent to the low frequency primary signal. 

2. The sea bottom does not affect the secondary field. This 
assumption is valid prov: ded that the sea depth d 2 is greater 
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than twice the low frequency skin depth of the sea, i.e., &2 > 2 
62 ( 1 °) • These two assumptions allow a halfspace inversion. The 
algebra and computer algorithms for the halfspace inversion are 
given in the Appendix. 

It is assumed that the altitude h is independently 
determined by a radar or Laser altimeter. Thus the inversion 
results in a local value for the sea ice thickness d^ and the 
conductivity of the sea water, a 2 ■ These results are employed in 
the inversion of the high frequency data to determine the sea ice 
conductivity . 

Outline of high frequency inversion 

First a halfspace inversion of the high frequency data 
is performed. This produces an effective skin depth 6 e ff which 
lies in the range 6 2 < (, eff — ^ 1 ’ an< ^ a function of ice 
thickness dj_. The effec ive high frequency skin depth is 
combined with the altimeter reading h to form the ratio A e ff = 
2h/6 e f f . The ad hoc normalization function employed in M1M 
inversion is a function of A e ££, i.e., 

MIM field - Normalized field = F(A eff ) Sommerfeld field . 
It turns out that the same normalization function F(A e ff) ma y he 
employed for both coil configurations. 

For d^ « 6 ^, 6 e f f ~ ^2 for d^ > 25^, then 6 ggg — 

5-j^ Since d]^ is known from the low frequency inversion, this 
latter case may be recognized and hence the first layer 
conductivity is determined from & e ff by 
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O I ” 2 /[/i 0 f ( hi) s eff 2 ] 

where /zq is the vacuum magnetic permittivity. The condition 
di > 2 61 only occurs for thick (d-^ > 10 m), highly conducting ( o ^ 

> 0.027 S/m) sea ice. For the intermediate, more general 
situation where 5 e ff/ d ^ 1 s of the order of unity, the inversion 
procedure to be used is that described below. 

The MIM relat ionship between the complex two-layer 
correction factor Q and the high frequency AEM field is 
algebraically transformed into two simultaneous transcendental 
real equations with argument d^/5^, where is the unknown 
quantity. All other quantities in these equations are known. 

Each of these equations has in general several roots, BUT only one 
common root. The explicit functions that occur respectively in 
these equations are tan(d^/ 6 |) and tarvh(d^/5 ^) . A root-finding 
algorithm is first applied to the tan(d^/ 6 ^) equation. When a 
root is determined, that root is inserted in the tanh(d^/ 6 ^) 
equation to test if it is also a root of the tanh(d^/ 6 -^> 
equation. If not, the algorithm continues in its determination 
of the real roots of the tan(d^/<5^) equation until the root is 
found that simultaneously satisfies both equations. The first 
layer skin depth 5 ^ and conductivity o ^ are given by that 
simultaneous root . 

The range of applicability of the root finding 
algorithm is given by 0.02 < dj/ 6 ]_ <2.5. These limits can be 
understood in physical terms. For d-^/5^ > 2.5 the sea ice is 
effectively a halfspace as has been already noted, and a two 
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layer model is inappropriate. For dy/6y < 0.02 the perturbation 
produced on the secondary AEM field by the sea ice cover is lost 
in the computer "noise" caused by roundoff, etc., and will 
certainly be undetectable in the noise and drift present in even 
ideal real data, where noise and drift are greater than about 1 
ppm. 

The lower ice thickness limit on the detectability of 
sea ice conductivity is illustrated in the following table which 
assumes a value of sea water conductivity of a 2 - 2.7 S/m and an 
operating frequency of 230 Khz. 


o 2 /oi 

5 1 

minimum d 

100 

-6m 

-0.1m 

1000 

- 20 m 

i 

o 

B 


The algebraic details of this procedure and the root finding 
algorithm are given in the Appendix. 

Results 

The MIM inversion procedure that has been described is 
applied to several sea ice models. In all of the models used the 
low and high frequencies assumed for the AEM system are 1 kHz and 
250 kHz, respectively; the altitude of the AEM bird is 25 m; the 
conductivity of the sea water a 2 is 2.7 S/m; and the conductivity 
of the sea ice ay for each model has input values of 0.027 S/m, 
0.0054 S/m, and 0.0027 S/m. Thus the ratio K of the 
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conductivities of sea wa er to sea ice has the values 100, 500, 
and 1000, respectively. 

With these general conditions, the first model of ice 
thickness versus range (fiducial number) is given in Figure 1. 

The ice thickness increases linearly with increasing range. 

Horizontal Coplanar Case 

The results of the inversion of the ZZ field for 
are shown in Figure 2. Ihe inversion values for o ^ are in fair 
agreement with the input values except for the case with a ^ — 
0.027 S/m. The problem occurs at an ice thickness of 
approximately 9.5 m. For = 0.027 S/m the skin depth of the 
sea ice is about 6 m, thus the ratio of ice thickness to skin 
depth (which is the argument of both the tan and tanh functions) 
is about pi/2, where the tangent becomes singular and double 
valued. More importantly, in the immediate vicinity of pi/2, 
tan(d^/6^) varies rapidly. In spite of this, the root finding 
inversion algorithm stil 1 works when the exact forward MIM field 
ZZ(MIM) is used as the input field (see Table 1). When a 
simultaneous root cannot be found for the normalized Sommerfeld 
field in the vicinity of tt/ 2, a value of 1.55 is assumed for x. 
See Table 2. It is the residual differences between the 
normalized Sommerfeld f eld (or real field data) and the exact 
MIM field that causes the root finding algorithm that we are 

currently using to fail for x =* ^\/&\ “ 7T /2 • 

It should be noted that this value of tt/2 will most 
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likely not be encountered in field surveys where ice 
conductivities will generally be less than 0.0054 S/m (K - 500). 
Table 3 shows that for K - 500, x is less than n/2 for sea ice 
thicknesses up to 20 m. 

A shallow ice keel model is shown in Figure 3. Figure 4 
shows the values of ct x lor this model produced by the inversion 
algorithm for K - 100, 500, and 1000. These results are also 

listed in Tables 5, 6, and 7. 

In all of the tables we have included the results of 

the half space inversion of the high frequency data which gives 

°e£f- It: can be S6en f ° r the C3Se K “ 100 When X > 2 ' 4 ’ ° eff = 

q , This demonstrates that when the ice thickness is greater 

input * 

than 2.4 skin depths, a halfspace inversion yields good results 
for the ice conductivity. Although values of x > 2.4 will 
probably not be found in survey data taken at a high frequency of 
250 kHz, still higher frequencies of about 1 MHz will bring x 
into this range. 

Finally, the results of the inversion are shown in 
Figure 3 for a shallow ice keel model. The tabulated results are 

shown in Tables 5, 6, and 7. 

In summary, the present inversion algorithm for a ^ 

works well except in the vicinity of d 1 /6 1 - V 2 - We are 
continuing efforts to modify and improve the existing algorithm. 

Vertical Coaxial Case 

Results similar to the horizontal coplanar case are 
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obtained for this coil configuration. 

Appendix 

The MIM representation of the normalized secondary 
field produced by induced ohmic currents in a two- layered 
conducting model (see Figure 1) for a horizontal coplanar coil 
pair is given by 

(H s /H p ) - ZZ - [2R 2 - 1 ] / [R 2 + 1] 5/2 (A- la) 

and for a vertical coaxial coil pair is given by 

(H s /H p ) - XX - [R 2 - 2]/[R 2 + 1] 5/2 , (A- lb) 

where for both cases 

R = [ 2h + (1-i) Q 6^/p , (A-2) 

and [2h + (1-i) Q i f: the complex vertical distance 

separating the primary dipole source from the image source, h is 
the real altitude of the bird above the first layer surface, 
[(1-i) Q 6]J/2 is the complex distance below the first layer 
surface of the image plane, p is the coil spacing, is the 
first layer (sea ice) skin depth, and finally Q is the two- layer 
correction factor given by 

Q . [{S x /s 2 ) + tanhl(1+i)di/«i)]/ll + («l/ 6 2> tanh{ (l+Ddj/fi].) ] 
d^ is the first layer thickness and S 2 is the skin depth of the 
second layer (sea water). The thickness or depth of the second 
layer, d 2 , is assumed to be greater than 2 S 2 in this analysis. 
(In order to determine the sea depth, a third lower frequency 
signal must be employed.) The R function can be interpreted 
geometrically as cotan j> , where <f> is the complex angle indicated 


(A- 3) 
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in Figure 1 . 


In MIM inversion schemes for both cases, Equation (A-2) 
is inverted by means of a polynomial expression which gives R as 
a function of ZZ or XX: 

l/R - as + bs 3 + cs 5 + ds 7 + es 9 , (A-4) 

where s - (ZZ/2)^/ 3 for the horizontal coplanar configuration and 
(2XX)^/ 3 for the vertical coaxial case. The values for the 
coefficients for the twc coil configurations are given in the 
following table. These Lnverted relationships are of paramount 
importance in all MIM inversion routines, i.e., the values of the 
model parameters are determined from the value of R calculated 
from the AEM fields. 


coil configuration 

a 

b 

c 

d e 

horizontal coplanar 

1 

1 

2.069 

-3 125.1 

vertical coaxial 

1 

1.5 

1.125 

32.08 203.7 


Half space inversions 



A halfspace inversion is defined by the condition Q 
1. Thus Equation (A-2) can be inverted to give values for the 
bird altitude h and the skin depth 6^ 

2h / p — R^ + 1*2 

S\/p ” ' r 2 

where Rj_ and R 2 are the real and imaginary components of R. In 
the event that the conditions which make Q - 1 are not satisfied 
(Q - 1 if o 2 - o x or d-, >26^, then the halfspace inversion 
results in an effective skin depth and bird altitude. 
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Two layer inversion for sea ice conductivity 


In this inversion it is assumed that the only unknown 
is the first layer skin depth, i.e., the altitude, the second 
layer skin depth, and the first layer thickness are known from 
the altimeter reading and the low frequency inversion results. 

The algebraic manipulations of Equations (A-2) and (A-3) that 
result in the two simultaneous equations mentioned in the 
narrative proceed as follows: 

First Q6^ is solved for explicitly from Equation (A-2) 

to produce 

Q6 1 - [p R(hi) - 2h]/(l - i) . (A-5) 

Next the expression for Q from Equation (A-3) is substituted into 
Equation (A-5) and the resultant equation is rearranged to solve 
explicitly for tanh{ (1+ i )d^/6^ ) : 

tanh{ ( 1+i) d^/S^ ) = [ (C • 1) (§ 2 /^ 1 ^ 3/tl * GS 2 /f>\ 

^ D(6 1 ,5 2 C) ■ (A-6) 

C is a known complex nunber given by 

C - [p R(hi) - 2 h]/I(l - i) « 2 1 
and D is a complex function of the unknown 6^ and the known 
quantities C and S 2 - If we set D - Dj_ + i D 2 and expand 
tanh{ (l+i)d^/5^ ) by 

tanh{ (l+i)d 1 /6 1 ) = [tanh(d 1 /« 1 ) + i tan(d 1 /6 1 ) ]/[ 1 + 

i tanh(dj_/5]_) tan(d^/6^)] , 

we get 

[tanh(d 1 /S 1 ) + i tan(d 1 /fi 1 )]/[l + i tanh(d 1 /5j) tan(d 1 /fi 1 )] 
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Dp + i D 2 


(A-7) 


Finally, if we equate the real and imaginary parts of Equation 
(A-7), we find 

tanh(dp/Sp) - Dp - D 2 tan(dp/6p) tanh(dp/6p) 
tan(dp/*p) “ D 2 + Dp tan(dp/6p) tanh(dp/Sp) 

These can be combined to give 

D 2 tan 2 (dp/Sp) - (Dp 2 + D 2 2 ' D tan(dp/Sp) ' D 2 “ 0 (A ' 8a) 

and 

Dl tanh 2 (dp/Sp) - (Dp 2 4 D 2 2 + D tanh(dp/«p) + Dp = 0 . (A-8b) 

Both of these equations are of the quadratic form ax + b x + c 
and hence explicit expressions for tanh(dp/«p) and tan(dp/«p) can 
be written. But first it should be noted that for the tan 
equation, a — - c - D 2 , and thus 

tan(dp/6p) - - b/2a ± t(b/2a) 2 + 1 1 1/2 . 

and similarly for the tanh equation, since a - c - D lt 
tanh(dp/6p) = - b/2a + [ (b/2a) - 1] 

The physical constraint that the roots of these 
equations be real and positive results in the following 


solutions : 

tan(dp/6p) =° Fp(dp/6p) * [Fp (dp/6p) + f) 
for 0 < (dp/«p) < Jf/2, and 
tan(dp/6p) - Fp - [Fp^ + 1] ^ ’ 

for w/2 < (dp/ 6 p) < if, where 

Fp = [Dp 2 + D 2 2 - l]/(2 D 2 ) 

tanh(dp/6p) = F 2 + [F 2 - 1] 


(A-9a) 


(A-9b) 


(A-10) 
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where 


F 2 = [D^ 2 4 D 2 2 + l]/(2 D-^ 

A root finding algorithm given below is used to find 
the real positive roots of Equation (A-9a) or (A-9b) . These 
roots are substituted into Equation (A-10) to find the one root 
of Equation (A-9a) or (A 9b) that is simultaneously a root of 
Equation (A-10). This value of (d is used to determine S 1 

and in turn 
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TABLE 1 


SIGMA1 (ACTUAL) =0 . 027 S/M 
WORKING WITH 


ZZMIM 

Sw R Sl ANi 5 s?G2?°SOLVE FOR X=D1/DELT1 AND HENCE SXG1 . 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVIT , 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 

X 

0.0 

0 . 1 

0.0156 

1.0 

0 . 5 

0.0781 

2.0 

1.0 

0.1563 

3 . 0 

1.5 

0.2344 

4.0 

2 . 0 

0.3250 

5.0 

2 . 5 

0.4063 

6.0 

3 . 0 

0.4875 

7 . 0 

3 . 5 

0.5687 

8.0 

4.0 

0.6503 

9.0 

4 . 5 

0.7312 

10 . 0 

5.0 

0.8125 

11.0 

5.5 

0.9000 

12 . 0 

6.0 

0.9818 

13 . 0 

6 . 5 

1.0636 

14 . 0 

7 . 0 

1.14 55 

15.0 

7 . 5 

1.2273 

16.0 

8 . 0 

1.309 1 

17 . 0 

8 . 5 

1 .3909 

18 . 0 

9 . 0 

1.4727 

19.0 

9 . 5 

1.5500 

20.0 

10.0 

1.6355 

21.0 

10 . 5 

1.7172 

22 . 0 

11.0 

1 .7990 

23 . 0 

11.5 

1.8808 

24.0 

12 . 0 

1.9626 

25.0 

12 . 5 

2 . 0443 

26.0 

13 . 0 

2 . 1218 

27.0 

13 . 5 

2.2035 

28.0 

14.0 

2.2851 

29 . 0 

14 . 5 

2 . 3667 

30.0 

15 . 0 

2.4433 

31.0 

15 . 5 

2 . 5239 

32.0 

16 . 0 

2.6115 

33 . 0 

16.5 

2.6931 

34 . 0 

17 . 0 

2.7747 

35.0 

17 . 5 

2 . 8563 

36.0 

18 . 0 

2.9379 

37.0 

18 . 5 

3 . 0195 

38.0 

19.0 

3 . 1012 

39 . 0 

19.5 

3 . 1828 

40.0 

20 . 0 

3 .2644 


SIG1 

SIGeff 

0.0247 

2 . 6960 

0.0247 

2.5941 

0 . 0247 

2.2813 

0.0247 

1.7051 

0 . 0268 

1.1303 

0.0268 

0.7022 

0.0268 

0 .4009 

0.0268 

0.2198 

0 . 0268 

0.1332 

0.0268 

0.0884 

0.0268 

0 . 0613 

0.0271 

0.0461 

0.0271 

0.0373 

0.0271 

0.0315 

0 .0271 

0 . 0277 

0.0271 

0.0253 

0.0271 

0.0238 

0.0271 

0.0228 

0 . 0271 

0.0224 

0.0270 

0 . 0221 

0.0271 

0.0222 

0.0271 

0.0223 

0 . 0271 

0.0226 

0.0271 

0.0230 

0.0271 

0.0234 

0.0271 

0.0237 

0.0270 

0.0241 

0.0270 

0.0244 

0.0270 

0 . 0247 

0.0270 

0.0250 

0 . 0270 

0 . 0252 

0.0270 

0.0254 

0.0270 

0.0255 

0.0270 

0 . 0257 

0 . 0270 

0.0258 

0.0270 

0.0258 

0.0270 

0.0259 

0.0270 

0.0259 

0.0270 

0.0259 

0.0270 

0.0259 

0.0270 

0.0259 




TABLE 2 


SIGMA1 (ACTUAL) =0.027 S/M 

WORKING WITH 

ZZSOM 


MD 5 S°IG2?°S0LV E FOR X=D1/DELT1 and hence SIG1. 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVIT , 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID D1 X 

0.0 0.1 0 . 0156 

1.0 0.5 0.0781 

2.0 1.0 0.1563 

3.0 1.5 0.2344 

4.0 2.0 0. 3375 

5.0 2.5 0.4087 

6.0 3.0 0.4904 

7.0 3.5 0.5722 

8.0 4.0 0 . 6664 

9.0 4.5 0.7497 

10.0 5.0 0.8330 

11.0 5.5 0.9163 

12.0 6.0 0.9996 

13.0 6.5 1.0329 

14.0 7.0 1.1662 

15.0 7.5 1.2495 

16.0 8.0 1.5500 

17.0 8.5 1.5500 

18.0 9.0 1.5500 

19.0 9.5 1. 4875 

20.0 10.0 1.5500 

21.0 10.5 1.7131 

22.0 11.0 1.7947 

23.0 11.5 1.8763 

24.0 12.0 1.9579 

25.0 12.5 2.0394 

26.0 13.0 2.1210 

27.0 13.5 2.2026 

28.0 14.0 2.2842 

29.0 14.5 2.3657 

30.0 15.0 2.4473 

31.0 15.5 2 . 5289 

32.0 16.0 2.6105 

33.0 16.5 2.6921 

34.0 17.0 2.7736 

35.0 17.5 2.8552 

36.0 18.0 2.9383 

37.0 18.5 3.0199 

38.0 19.0 3.1024 

39.0 19.5 1.5500 

40.0 20.0 1.5500 


SIG1 

SIGeff 

0.0247 

2.6973 

0.0247 

2.5953 

0 . 0247 

2 . 2816 

0 . 0247 

1.7033 

0.0289 

1.1265 

0 . 0271 

0.6976 

0 . 0271 

0.3973 

0 . 0271 

0.2178 

0.0281 

0.1324 

0 . 0281 

0.0884 

0.0281 

0 . 0619 

0.0281 

0.0471 

0 . 0281 

0.0386 

0.0281 

0.0330 

0 .0281 

0.0294 

0.0281 

0.0270 

0.0380 

0.0256 

0.0337 

0 . 0247 

0.0301 

0.0242 

0. 0248 

0.0240 

0.0243 

0.0240 

0 . 0270 

0.0241 

0.0270 

0.0244 

0.0270 

0.0247 

0.0270 

0.0250 

0.0270 

0.0253 

0.0270 

0.0256 

0.0270 

0.0259 

0.0270 

0.0261 

0 . 0270 

0 . 0264 

0 . 0270 

0.0265 

0 . 0270 

0.0267 

0 . 0270 

0.0268 

0.0270 

0 . 0269 

0 .0270 

0.0270 

0.0270 

0.0271 

0 . 0270 

0.0271 

0 . 0270 

0.0271 

0.0270 

0.0271 

0.0064 

0.0271 

0 . 0061 

0.0271 




ABLE 3 


SIGMA1 (ACTUAL) =0 . 0054 S/M 

WORKING WITH 

ZZSOM 

XSTART= 0.5000000 

KJIOW D1 AND SIG2, SOLVE FOR X=D1/DELT1 AND HENCE SIG1 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVITY, 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 

X 

SIG1 


SIGeff 

0.0 

0.1 

0.0078 

0.0062 


2 . 6998 

1.0 

0.5 

0.0391 

0.0062 


2 . 6784 

2.0 

1.0 

0.0781 

0. 0062 


2.6794 

3.0 

1 . 5 

0. 1172 

0 . 0062 


2 . 5616 

4.0 

2.0 

0.1563 

0.0062 


2.3817 

5.0 

2.5 

0.1953 

0.0062 


2 . 1671 

6.0 

3 . 0 

0.2344 

0 . 0062 


1.7031 

7.0 

3 . 5 

0.2734 

0 . 0062 


1.1052 

8.0 

4.0 

0.2930 

0 . 0054 


0.7481 

9.0 

4 . 5 

0.3296 

0.0054 


0.5329 

10.0 

5.0 

0.3662 

0.0054 


0.3495 

11.0 

5.5 

0.4028 

0.0054 


0.2399 

12.0 

6.0 

0.4395 

0.0054 


0.1807 

13.0 

6.5 

0.4886 

0.0057 


0.1303 

14 . 0 

7.0 

0. 5262 

0.0057 


0.0979 

15.0 

7 . 5 

0.5637 

0.0057 


0.0723 

16.0 

8 . 0 

0.6076 

0 . 0058 


0.0536 

17.0 

8 . 5 

0.6455 

0.0058 


0.0414 

18.0 

9.0 

0.6835 

0.0058 


0.0325 

19.0 

9.5 

0.7215 

0.0058 


0 . 0264 

20.0 

10.0 

0.7595 

0 . 0058 


0.0213 

21.0 

10.5 

0.7974 

0 . 0058 


0.0178 

22.0 

11.0. 

0.8354 

0.0058 


0.0148 

23.0 

11.5 

0.8734 

0.0058 


0.0126 

24.0 

12.0 

0.9114 

0.0058 


0.0109 

25.0 

12 . 5 

0.9493 

0.0058 


0.0097 

26.0 

13 . 0 

0.9873 

0.0058 


0.0087 

27.0 

13 . 5 

1.1253 

0.0070 


0.0082 

28.0 

14.0 

1.1670 

0.0070 


0.0075 

29.0 

14 . 5 

1.2086 

0 . 0070 


0.0070 

30.0 

15.0 

1 . 2503 

0.0070 


0.0066 

31.0 

15.5 

1.2112 

0.0062 


0.0063 

32.0 

16.0 

1.2503 

0.0062 


0.0060 

33.0 

16.5 

1.2894 

0.0062 


0.0058 

34.0 

17.0 

1.3285 

0.0062 


0.0056 

35.0 

17.5 

1.3675 

0.0062 


0.0055 

36.0 

18 . 0 

1.4066 

0.0062 


0.0054 

37.0 

18.5 

1 . 5500 

0 . 0071 


0.0053 

38.0 

19.0 

1.5500 

0.0067 


0.0052 

39.0 

19.5 

1.5500 

0.0064 


0.0052 

40.0 

20.0 

1.5500 

0.0061 


0.0051 




TABLE 4 


SIGMA1 (ACTUAL) -0 . 0027 S/M 
WORKING WITH 


ZZSOM 
XSTART= 
KNOW D1 
SIGef f 
AND IT 


0 . 5000000 
AND SIG2 , S 
IS THE HALF- 
IS COMPUTED 


OLVE FOR X=D1/DELT L AND HENCE SIG1. 
SPACE EFFECTIVE CONDUCTIVITY, 

ONLY WHEN INVERTING ZZSOM 


FID 

0.0 

1.0 

2.0 

3.0 

4.0 

5.0 

6.0 

7 . 0 

8 . 0 
9.0 

10.0 
11.0 
12 . 0 

13 . 0 

14 . 0 

15.0 

16.0 

17 . 0 

18. 0 

19 . 0 

20.0 
21.0 
22 . 0 

23 . 0 

24 . 0 

25.0 

26.0 

27 . 0 

28 . 0 

29.0 

30 . 0 

31.0 

32 . 0 

33 . 0 

34 . 0 

35.0 

36 . 0 

37 . 0 

38 . 0 

39 . 0 

40.0 


D1 X 

0.1 0.0078 

0.5 0.0391 

1.0 0.0781 

1.5 0.0586 

2.0 0.0781 

2.5 0.0977 

3.0 0.1672 

3.5 0.1951 

4.0 0.2229 

4.5 0.2508 

5.0 0.2786 

5.5 0.2871 

6.0 0.3135 

6.5 0.3395 

7.0 0.3657 

7.5 0.3913 

8.0 0 . 4180 

8.5 0.4566 

9.0 0.4835 

9.5 0.5103 

10.0 0 . 5372 

10.5 0.5640 

11.0 0.5909 

11.5 0.6177 

12.0 0.6446 

12.5 0.6715 

13.0 0.6983 

13.5 0.7252 

14.0 0.7520 

14.5 0.7789 

15.0 0.8058 

15.5 0.9326 

16.0 0.9627 

16.5 0.9928 

17.0 1.0229 

17.5 1.0529 

18.0 1.0830 

18.5 1.1131 

19.0 1.1432 

19.5 1.1733 

20.0 1.2034 


SIG1 
0 . 0062 
0.0062 
0 . 0062 
0 . 0015 
0.0015 
0.0015 
0 . 0031 
0 . 0031 
0 . 0031 
0 . 0031 
0 .0031 
0 . 0028 
0 . 0028 
0 . 0028 
0 . 0028 
0 . 0028 
0.0028 • 
0 . 0029 
0 . 0029 
0 . 0029 
0.0029 
0 . 0029 
0 . 0029 
0 . 0029 
0 . 0029 
0 . 0029 
0 . 0029 
0 . 0029 
0 . 0029 
0.0029 
0 . 0029 
0 . 0037 
0.0037 
0.0037 
0 . 0037 
0 . 0037 
0.0037 
0.0037 
0 . 0037 
0.0037 
0.0037 


SIGeff 
2.7001 
2.6894 
2.7387 
2.7221 
2.7030 
2.7093 
2 . 3818 
1.6995 
1.2689 
1.0030 
0.7038 
0.5150 
0.4217 
0.3189 
0.2518 
0.1897 
0.1405 
0.1086 
0.0843 
0.0680 
0.0528 
0 . 0426 
0.0332 
0.0265 
0.0217 
0.0180 
0.0151 
0 . 0128 
0.0109 
0 . 0095 
0.0083 
0 . 0073 
0.0067 
0 . 0061 
0 . 0055 
0 . 0051 
0.0047 
0.0044 
0.0041 
0 . 0038 
0.0036 




TABLE 5 


SIGMA1 (ACTUAL) =0.027 S/M 

WORKING WITH 

ZZSOM 


$flOW R Dl AND 5 SIG2?°S0LVE FOR X=D1/DELT1 AND HENCE SIG1 
SIGeff IS THE HAI.F-SPACE EFFECTIVE CONDUCTIVIT , 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 

X 

0.0 

0.1 

0.0156 

1.0 

0.5 

0.0781 

2.0 

1.0 

0.1563 

3.0 

1.5 

0.2469 

4 . 0 

2 . 0 

0.3292 

5.0 

2.5 

0.4115 

6.0 

3 . 0 

0.4938 

7 . 0 

3 . 5 

0 . 5760 

8 . 0 

4.0 

0.6533 

9.0 

4.5 

0.7469 

10 . 0 

5.0 

0.8299 

11 . 0 

5.5 

0.9191 

12 . 0 

6.0 

1.0027 

13 . 0 

6.5 

1.0862 

14 . 0 

7 . 0 

1.1652 

15 . 0 

7 . 5 

1.2484 

16.0 

8 . 0 

1.5 5 00 

17.0 

8 . 5 

1.5500 

18 . 0 

9 . 0 

1.5500 

19.0 

9.5 

1.5000 

20.0 

10.0 

1.6274 

21.0 

9.5 

1.5500 

22 . 0 

9 . 0 

1. 5500 

23 . 0 

8 . 5 

1 .5500 

24 . 0 

8 . 0 

1.5500 

25.0 

7.5 

1.2500 

26.0 

7.0 

1. 1667 

27.0 

6.5 

1.0833 

28.0 

6.0 

1.0000 

29 . 0 

5.5 

0.9L67 

30.0 

5 . 0 

0.8333 

31.0 

4 . 5 

0.7500 

32 . 0 

4 . 0 

0 . 6667 

33.0 

3.5 

0.5833 

34 . 0 

3 . 0 

0 . 5000 

35.0 

2 . 5 

0 .4167 

36.0 

2 . 0 

0.3333 

37.0 

1.5 

0 .2500 

38 . 0 

1.0 

0. 1667 

39.0 

0.5 

0.0333 

40.0 

0.1 

0. 0167 


SIG1 

SIGeff 

0.0247 

2 . 6973 

0.0247 

2 . 5953 

0 . 0247 

2.2046 

0 . 0274 

1.6041 

0.0274 

1.0232 

0.0274 

0.6057 

0 . 0274 

0.3527 

0.0274 

0 .2107 

0.0274 

0.1324 

0 . 0279 

0.0884 

0.0279 

0.0631 

0 . 0283 

0. 0481 

0 . 0283 

0.0388 

0.0283 

0.0330 

0.0281 

0 . 0294 

0.0281 - 

0.0270 

0.0380 

0.0256 

0. 0337 

0.0247 

0.0301 

0 . 0242 

0 . 0253 

0.0240 

0.0268 

0.0241 

0.0270 

0.0240 

0.0301 

0.0242 

0 . 0337 

0.0247 

0.0380 

0.0256 

0.0281 

0.0270 

0.0281 

0.0294 

0.0281 

0.0330 

0.0281 

0.0388 

0 . 0281 

0.0481 

0 . 0281 

0.0631 

0.0281 

0.0884 

0.0281 

0.1324 

0.0281 

0.2107 

0.0281 

0 .3527 

0.0281 

0.6057 

0.0281 

1.0232 

0 . 0281 

1.6041 

0.0281 

2.2046 

0.0281 

2.5953 

0.0281 

2.6973 




TABLE 6 


SIGMA1 (ACTUAL) =0.0054 S/M 
WORKING WITH 


Z Z SOM 

A^S^r SOLVE FOR X-D1/DELT1 AND HENCE SIG1. 
SIGeff IS THE HALF-SPACE EFFECTIVE CONDUCTIVITY, 

AND IT IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 

X 

0.0 

0.1 

0.0078 

1.0 

0.5 

0.0391 

2 . 0 

1.0 

0.0781 

3.0 

1.5 

0.1172 

4 . 0 

2 . 0 

0.1563 

5.0 

2.5 

0.1953 

6.0 

3 . 0 

0.2344 

7.0 

3.5 

0.2734 

8 . 0 

4 . 0 

0 . 2930 

9.0 

4.5 

0.3296 

10.0 

5.0 

0.3662 

11.0 

5.5 

0 .4026 

12 . 0 

6.0 

0.4395 

13 . 0 

6.5 

0.4886 

14 . 0 

7 . 0 

0 .52 62 

15 . 0 

7 . 5 

0.5637 

16.0 

8 . 0 

0 . 6076 

17 . 0 

8.5 

0. 645.5 

18 . 0 

9 . 0 

0.68 3 5 

19 . 0 

9.5 

0.7215 

20.0 

10.0 

0.7 59 5 

21.0 

9.5 

0.7215 

22.0 

9.0 

0.6835 

23 . 0 

8 . 5 

0.6500 

24 . 0 

8 . 0 

0.6113 

25.0 

7 . 5 

0.5646 

26.0 

7 . 0 

0.5269 

27 . 0 

6 . 5 

0.4893 

28.0 

6 . 0 

0.4517 

29 . 0 

5.5 

0.4140 

30 . 0 

5 . 0 

0.3764 

31.0 

4 . 5 

0.3387 

32 . 0 

4 . 0 

0.3011 

33.0 

3 . 5 

0.2635 

34.0 

3 . 0 

0.2258 

35.0 

2 . 5 

0.1862 

36.0 

2 . 0 

0.1506 

37.0 

1.5 

0.1129 

38 . 0 

1.0 

0.07 5 3 

39.0 

0.5 

0.0 37 6 

40.0 

0.1 

0.00^5 


SIG1 

SIGeff 

0 . 0062 

2 . 6998 

0.0062 

2 . 6784 

0.0062 

2 . 5848 

0 . 0062 

2 . 3914 

0.0062 

2 . 0990 

0 . 0062 

1.7409 

0.0062 

1.3686 

0 . 0062 

1.0287 

0.0054 

0.7481 

0 . 0054 

0.5329 

0 . 0054 

0.3759 

0 . 0054 

0.2650 

0 . 0054 

0. 1879 

0.0057 

0.1348 

0 . 0057 

0.0979 

0 . 0057 

0 . 0723 

0 . 0058 

0 . 0544 

0.0058 

0.0417 

0.0058 

0.0325 

0 . 0058 

0.0259 

0 . 0058 

0 . 0210 

0.0058 

0.0259 

0 . 0058 

0.0325 

0.0059 

0.0417 

0.0059 

0.0544 

0.0057 

0.0723 

0.0057 

0.0979 

0.0057 

0.1348 

0 . 0057 

0.1879 

0.0057 

0.2650 

0.0057 

0.3759 

0 . 0057 

0.5329 

0 . 0057 

0.7481 

0 . 0057 

1.0287 

0.0057 

1.3686 

0 . 0057 

1.7409 

0.0057 

2.0990 

0.0057 

2.3914 

0.0057 

2.5848 

0.0057 

2.6784 

0 . 0057 

2 . 6998 




TABLE 7 


SIGMA1 (ACTUAL) =0.0027 S/M 
WORKING WITH 


ZZSOM 
XSTART= 
KNOW D1 
SIGef f 
AND IT 


Am 5 s?G2?°EOLVE FOR X-D1/DELT1 AND HENCE SKI. 
IS THE HALF-SPACE EFFECTIVE CONDUCTIVITY, 

IS COMPUTED ONLY WHEN INVERTING ZZSOM 


FID 

D1 X 

0.0 

0.1 0.0078 

1.0 

0.5 0.0391 

2.0 

1.0 0.0781 

3 . 0 

1.5 0 . 0586 

4.0 

2.0 0.0781 

5.0 

2.5 0.0977 

6.0 

3.0 0.1672 

7.0 

3.5 0.1951 

8 . 0 

4.0 0.2229 

9.0 

4.5 0.2508 

10.0 

5.0 0.2786 

11.0 

5.5 0.2874 

12 . 0 

6.0 0.3135 

13.0 

6.5 0.3396 

14 . 0 

7.0 0.3657 

15.0 

7.5 0.3918 

16.0 

8.0 0.4180 

17 . 0 

8.5 0.4566 

18.0 

9.0 0.4835 

19 . 0 

9.5 0.5103 

20.0 

10.0 0.5372 

21.0 

9.5 0.5103 

22 . 0 

9.0 0.4835 

23 . 0 

8.5 0 . 4 5 6 <3 

24 . 0 

8.0 0.4297 

25.0 

7.5 0.4029 

26.0 

7.0 0.3760 

27.0 

6.5 0.3492 

28 . 0 

6.0 0.3223 

29.0 

5.5 0.2954 

30.0 

5.0 0.2686 

31.0 

4.5 0.2417 

32.0 

4.0 0.2149 

33.0 

3.5 0.1880 

34 . 0 

3.0 0. 1612 

35.0 

2.5 0.1343 

36.0 

2.0 0 . 107 4 

37.0 

1.5 0.0806 

38 . 0 

1.0 0.0537 

39.0 

0.5 0.0269 

40.0 

0.1 0.0054 


SIG1 

SIGeff 

0 . 0062 

2.7001 

0 . 0062 

2 . 6894 

0.0062 

2 . 6414 

0.0015 

2 . 5379 

0.0015 

2 . 3690 

0.0015 

2.1376 

0.0031 

1.8600 

0 . 0031 

1 . 5616 

0.0031 

1 . 2689 

0.0031 

1 . 0030 

0 . 0031 

0.7760 

0 . 0028 

0 . 5912 

0 . 0028 

0 .4462 

0 .0028 

0.3353 

0.0028 

0 . 2518 

0 . 0028 

0.1897 

0.0028 

0.1437 

0.0029 

0.1096 

0 . 0029 

0 . 0843 

0.0029 

0 . 0654 

0 . 0029 

0.0513 

0.0029 

0.0654 

0.0029 

G .0843 

0 . 0029 

0.1096 

0.0029 

0.1437 

0.0029 

0.1897 

0.0029 

0.2518 

0 .0029 

0.3353 

0.0029 

0.4462 

0.0029 

0.5912 

0.0029 

0.7760 

0.0029 

1.0030 

0.0029 

1.2689 

0.0029 

1.5616 

0.0029 

1.8600 

0.0029 

2.1376 

0.0029 

2 .3690 

0 . 0029 

2.5379 

0.0029 

2.6414 

0.0029 

2 . 6894 

0 . 0029 

2.7001 




FIGURE 1 - Appendix 
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